Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm struggling with a small data set and a slow calculation. I have hundreds of small 2D data arrays and need to integrate across several lines parallel to the y-axis.

Let me start with the data array:

data={{51.6,10},{16.52,18},{48.835,12},{72.975,7},{61.085,19},{78.362,11},{24.312,25},
      {26.818,5},{37.336,5},{70.647,18},{48.783,35},{40.379,35},{75.926,8},{23.226,8},
      {54.727,18},{116.046,12},{20.597,20},{21.973,28},{74.821,7},{56.408,51}};

After rescaling and truncating the distribution, it looks like this:

dnorm=Transpose[Rescale/@Transpose[data]];
pdf[x_,y_]=PDF[TruncatedDistribution[{{0,∞},{0,∞}},
               SmoothKernelDistribution[dnorm,Automatic,{"Radial","Gaussian"}]],{x,y}];
ContourPlot[pdf[x,y],{x,0,1},{y,0,1}]

plot of distribution

I am interested in the profile across a line parallel to the y-axis (say e.g. $x=0.5$) and do handle this profile as a distribution for this specific x-value. I need the 95% quantile for this specific distribution at the specific x-value.

Finally, I am interested in all the quantile-values for all x-values between 0 and 1 (in the rescaled domain) and do rescale them back afterwards.

I have written the following function, which works but need approx. 10 seconds (even without the line plot and sampling at 20 data points between 0 and 1). I have to evaluate hundreds of data sets with this function, so I need to speed it up.

myquantile[d_List,q_]:=Block[{dnorm,pdf,x,y,func,int,ξ,ψ,sol,invx,invy,lp},
  dnorm=Transpose[Rescale/@Transpose[d]];
  pdf[x_,y_]=PDF[TruncatedDistribution[{{0,∞},{0,∞}},
                 SmoothKernelDistribution[dnorm,Automatic,{"Radial","Gaussian"}]],{x,y}];
  func[ξ_]:=Block[{},
    int=NIntegrate[pdf[ξ,y],{y,0,1}];
    ψ/.FindRoot[NIntegrate[pdf[ξ,y],{y,0,ψ}]/int==q,{ψ,0.8,0.9},Evaluated->False]];
    {invx[x_],invy[x_]}=Rescale[x,{0,1},#]&/@({Min[#],Max[#]}&/@Transpose[d]);
    sol=Interpolation[Table[{invx[x],invy[func[x]]},{x,0,1,0.05}]
  ];
  lp=Plot[sol[x],{x,invx[0],invx[1]},AxesOrigin->{invx[0],invy[0]}];
  {sol,lp}
];
myquantile[data, 0.95] // Timing

The final InterpolatingFunction looks like this (in the regular domain):

plot of interpolant

Can anybody point me in the right direction? Is it possible to get the line distribution in a different and faster way? A MarginalDistribution is pretty much what I want, but it is the distribution for all x-values and it gives me no information on individual x-values.

Thanks!

share|improve this question
1  
@0x4A4D: how have you changed the [Psi] etc. into the "real" symbols? It was showing up differently, when I posted this question? –  akm Jun 20 '13 at 13:21
    
Well, magic... okay, no. I knew the Unicode for these, so I put 'em in. –  J. M. Jun 20 '13 at 13:24

1 Answer 1

up vote 1 down vote accepted

You can setup a partial differential equation to integrate along $y$. It costs less than 0.1 second on my old machine.

sol = NDSolve[{
               D[cdf[x, y], y] == pdf[x, y],
               cdf[x, 0] == 0
              },
              cdf, {x, 0, 1}, {y, 0, 1}]

quantileLine = ContourPlot[Evaluate[
                      cdf[x, y]/cdf[x, 1] == .9 /. sol[[1]]],
                      {x, 0, 1}, {y, 0, 1}]

enter image description here

Or if you want an explicit InterpolatingFunction:

quantileLineFunc = Interpolation@
    Cases[quantileLine,
        GraphicsComplex[pts_, others__] :>
            Part[pts,
             Cases[others, Line[idx_] :> idx, ∞][[1]]
             ]
        ][[1]]
share|improve this answer
    
@ Silvia: thanks a lot; your answer is very obvious, since you just use the dependency between the cdf and pdf. I was completely blocked in mind and haven't thought about other directions to solve this. Again, thanks a lot! –  akm Jun 20 '13 at 15:31
    
@akm You're welcome! –  Silvia Jun 20 '13 at 15:43

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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