How to simulate annual variability in bimodal climate data?

Question

I have some radiation data with a distribution that looks somewhat bimodal. I am looking for a distribution to fit the data (in a negative log likelihood function which will be minimized using optim()). I have tried fitting a normal and a beta binomial without much success. I would prefer not to transform my data. Is there a distribution is flexible enough to fit with bimodal data?

Here is a histogram of my radiation data:

I have several climate variables: radiation, temperature and precipitation. I am creating a model for each of them so that I can simulate data from them. I have successfully fit a model to the temperature data using a normal distribution. I have also successfully fit a model to precipitation data using an exponential distribution. The only one I have not been able to fit is radiation.

I am fitting a sine curve to each data set because of the seasonal nature of it, and I am optimizing the curve using the negative log likelihood. The red curve is my rough estimate of where the fitted sine curve should be.

Here is all my radiation data plotted against Julian day on the x-axis:

Here is what I am ultimately trying to do: I want to simulate data at the high and low ends of the data cloud (i.e. years with lots of high temperatures and years with lots of low temperatures) to test the sensitivity of a different complex model I am using. I was able to use the normal distribution with temperature data and the exponential with precipitation, but I can't find one to work with radiation data given the slightly bimodal nature of the data.

Answer 1

Given that your goal is to seasonally detrend the data, your regression model will be something like $$R[t]=f[t]+\epsilon$$ where $R[t]$ is the radiation time series, $f$ is the (sinusoidal) trend, and $\epsilon$ is the residual (detrended) data.

Your question is about the (marginal) distribution $p[R]$, which appears to be bimodal. However this does not necessarily imply that the (conditional) distribution of the residuals $p[\epsilon]$ is bimodal!

In fact, if you plot a histogram for your sinusoidal trend $f[t]$, you should find that it is bimodal all by itself (with modes corresponding to the peak and trough).

I would suggest trying a simple sine curve fit* to the data, then look at the distribution of the residuals. (*For example, use whatever linear regression routine, but pass in $\sin[2\pi t/T]$ instead of $t$; assuming your period $T$ is known, e.g. 1 year.)

---

Update: Here is a simple example to illustrate the above ideas. (This should not be taken as a full solution!)

The figures below correspond to $$f[t] = 38\times\exp\left[0.3\times\sin\left(\tfrac{2\pi\,t}{365}\right)\right]-18$$ and $\epsilon = 2.5 z$ where $z$ is white noise with a standard normal distribution (total number of points is $n=1001$, on a uniform grid over $t\in[0,365]$).

This gives the following time series below

The corresponding (bimodal) histogram is

Notes:

• I put in the exponential simply to break the symmetry of the two modes.
• For this model the data $R[t]$ cannot simply be log-transformed.
• This example also assumes your data is censored (i.e. "$35$" really means "$\geq 35$").

Finally, the model form above could be fit to your data using nonlinear least squares. (Really there is 1 nonlinear parameter, the coefficient inside the exponential, so the least squares problem is "separable".)

Answer 2

My impression of your data (as provided) is of very weak bimodality. We don't have much context, but based on what is given, I'd say that attempting to define a bimodal distribution then fit it using some minimization routine would be less advisable than some data transformation...that is, if you even need one. The histogram certainly doesn't conform to normality but I've seen far worse. If you can provide more background we can be of more help.