We propose a method that allows for a rigorous statistical analysis of neural responses to natural stimuli, which are non-Gaussian and exhibit strong correlations. We have in mind a model in which neurons are selective for a small number of stimulus dimensions out of the high dimensional stimulus space, but within this subspace the responses can be arbitrarily nonlinear. Therefore we maximize the mutual information between the sequence of elicited neural responses and an ensemble of stimuli that has been projected on trial directions in the stimulus space. The procedure can be done iteratively by increasing the number of directions with respect to which information is maximized. Those directions that allow the recovery of all of the information between spikes and the full unprojected stimuli describe the relevant subspace. If the dimensionality of the relevant subspace indeed is much smaller than that of the overall stimulus space, it may become experimentally feasible to map out the neuron's input-output function even under fully natural stimulus conditions. This contrasts with methods based on correlations functions (reverse correlation, spike-triggered covariance,...) which all require simplified stimulus statistics if we are to use them rigorously.