0
$\begingroup$

I am trying to recreate some results from a simulation that I have but I am having trouble re-creating the expressions/funcitons in mathematica.

In my simulation I have a data set consisting of (z,x,p) points. I am interested in finding the the correlation between x and p, which I can easily compute by finding the average of their product . Below I have included plots of each the x, p and x*p data set as a function of z:

enter image description here enter image description here enter image description here

Now, I would like to analyze this system in a more theoretical/analytical way by expressing and/or fitting the data of the x and p columns as analytical expression/functions. The p function is a continuous function which I have been successful in describing. The x column on the other hand is what I am having difficulty with.

The x column can be easily visualized, but I am having a lot of trouble describing it mathematically in terms of an equation that I can multiply with p[z] and find the average. Describing the x function: The x function is gaussainly distributed along both the z-axis and x-axis (refer to above plot of x[z]). Meaning that most of the data points (or density if we're talking about a continuous function) are found at the core of the plot (0,0) and the density falls off as you move along z or x.

After, generating the x function I want to multiply it by the continuous p[z] and compute the average.

Any help would be greatly appreciated! I've been trying this for a while!

Some code (doesn't work) to help visualize my thought process:

sigmaX0 = 10^-4;(*define roughly the max sigmaX of the data at z=0*)
sigmaz = 50*10^-5;(*define roughly the sigma of the z axis*)

sigmax[z] := sigmaX0/(Sqrt[2 Pi] sigmaz) Exp[-(z^2/(2*(sigmaz)^2))](*define
the sigmax as a function of z; x is gaussianly dist across z*)


productOfXandP[z_] := NIntegrate[1/(Sqrt[2 Pi] sigmax[z]) Exp[-(x^2/(2*(sigmax[z])^2))]*p(z), {x, -Infinity, Infinity}]

average=NIntegrate[productOfXandP[z],{z,-Infinity,Infinity}]
$\endgroup$
  • 1
    $\begingroup$ You can use the build in multivariate normal distribution (see here) for the first part of your problem. $\endgroup$ – Armin Oct 4 '16 at 22:52
  • $\begingroup$ Maybe you could use RandomVariate with a NormalDistribution function in order to generate a table of normally distributed random values of x. However, it is not clear to me what is the correlation between x and z in this newly defined function of yours. Could you clarify that much? $\endgroup$ – FACamargo Oct 4 '16 at 23:14
  • 1
    $\begingroup$ that is probably just me being stupid... but what do you mean by average (along x?, along z?, along x- and z?)? If you would write down some formulas I might be of more use... $\endgroup$ – Armin Oct 4 '16 at 23:25
  • $\begingroup$ I edited my post to include some code to visualize my idea, and clarified the definition of average: <x*p>...essentially the average of the product of x and p. $\endgroup$ – user1886681 Oct 5 '16 at 0:31
  • 1
    $\begingroup$ there is an error in sigmax[z]. It needs to be sigmax[z_].... and there is also that second integral :) you should think about that one... it evaluates to p(z) by definition... ;) $\endgroup$ – Armin Oct 5 '16 at 0:59
2
$\begingroup$

I only try to give some suggestions here since I am far from being an expert in fitting or data analysis with mathematica (or any other software). Furthermore, at the end of that post, I'll give some topics which might be of interest for you.

@user1886681 since we had some problems in understanding each other I will summarize what I understood so far.

1.) You have two pdf's $p(z)$ and $X(x,z)$. You have an analytical expression for $p(z)$ and want to get one for $X(x,z)$ for further analysis.
2.) You guessed $X(x,z)$ is a bivariate normal distribution (see here) but you don't know how two obtain the moments $\mu_i$, $\sigma_i$ and the correlation.

The main idea is shortly summarized in the following commands

(* creating data for a distribution similar to the described one \
(bivariate normal distribution)*)
Xdata = RandomVariate[
   MultinormalDistribution[{0, 
     0}, {{0.25, 0.1*0.5*2}, {0.1*0.5*2, 4}}], 10000];

The created distribution $X(x,z)$ looks like:

(* show the distribution *)
Histogram3D[Xdata, 
 AxesLabel -> {"x", "y", "X(x,y)"}]

Similar to your distribution?

You want to get the parameters:

(* mean (Subscript[\[Mu], 1],Subscript[\[Mu], 2]),  \[Rho] is the \
correlation, Subscript[\[Sigma], i]=Sqrt[Subscript[variance, i]] *)
\
FindDistributionParameters[Xdata, 
 MultinormalDistribution[{mu1, 
   mu2}, {{sigmaX^2, rho*sigmaX*sigmaZ}, {rho*sigmaX*sigmaZ, 
    sigmaZ^2}}]]

Output:

{mu1 -> -0.0027094, mu2 -> -0.0103807, sigmaX -> -0.499168, 
 rho -> -0.0943205, sigmaZ -> 1.96191}  

There are some general suggestions for things you should clarify before further analyzing your data and your model:
1.) fit specific (other people in this forum will probably be able to help you)
1.1) How do I get the goodness of the fit?
1.2) Which method do I have to use?
1.3) How do I get the error of the estimated parameters?
1.4) Do I have to normalize my data $X(x,z)$? (I would guess so)

2.) problem specific (probably some other physicist)
2.1) what is a sensible distribution $X(x,z)$

3.) statistics (let $p(z)$ be some distribution)
3.1) You confused normalization and mean (average)
$\int_{-\infty}^{\infty} p(z) dz$ will give you some kind of normalization
$\int_{-\infty}^{\infty} z p(z) dz=\overline{z}$ will give you a mean
(that lead to some confusion regarding your definition of productOfXandP[z_])
3.2) correlation and how to calculate it (mathematica has a bulid in function for that)

$\endgroup$
  • $\begingroup$ Awesome thanks! I will look at this in detail... $\endgroup$ – user1886681 Oct 6 '16 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.