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 have the following expression, only slightly complicated, and I'd like to do two things:

  1. Plot the elements of the lists (though messy, all are in the same range) in terms of the range of the random variable in Expectation. (0.1 in the working example)
  2. Let the Plot work with Manipulate so that the user can experiment with different underlying values.

With all these SetDelayed used, I could not get Manipulate to work (Plot is simply slow). Basically, I am not sure if I can simply feed a new value for underlying into the last line (Expectation), or I need Replace, or some rule, a Block, or With?

underlying := {0.8,0.7,0.5,0.4,0.8,0.7}
answers := MapThread[(Min[1,Max[#1 + #2 ,0]])&,{underlying, {e1,e2,e3,e4,e5,e6}}]
disc :={answers[[5]],answers[[6]],answers[[3]]/answers[[1]],answers[[4]]/answers[[2]]} 
DB :={disc[[2]]/disc[[1]],disc[[1]]^2/disc[[2]]} 
db :={disc[[4]]/disc[[3]],disc[[3]]^2/disc[[4]]} 
deltasbetas = Expectation[{disc,DB,db},{e1,e2,e3,e4,e5,e6}
      \[Distributed]UniformDistribution[Table[{0,0.1},{6}]]]

(* {{0.85,0.75,0.647807,0.600891},{0.883373,0.965877},{0.931221,0.705072}} *)
share|improve this question
    
Could you post your (non-working) Plot? –  belisarius Jul 13 '12 at 17:38
    
@belisarius: Thanks for your code below, I tried to comment on it. –  László Jul 15 '12 at 18:17
add comment

3 Answers 3

Several observations as a basis for paper-and-pencil approach which, in turn, "might" suggest the first steps of an alternative formulation working with univariate expectations:

  • the components of the random vector e and, hence, those of the random vector answers are independent random variables
  • the remaining random variables discs,DB db are products of independent random variables, hence their expectations are the products of the expectations of constituent RVs
  • For each component of answers, we need four central moments (namely, 1,2, -1 and -2) to compute the expectations of the remaining random variables.

So, let

z[a_, e_] := Min[1, Max[a + e, 0]]
moment[a_, d_, m_Integer] := 
Expectation[z[a, e]^m, e \[Distributed] UniformDistribution[{0, d}]]

Table[{i, moment[a, d, i]}, {i, {1, 2, -1}}] // TableForm

enter image description here

EDIT: Using the observations above and a modification of Istvan's answer for making the legends:

expVal[a_?NumericQ, d_?NumericQ, m_?NumericQ] :=  
     NIntegrate[Min[1, Max[a + x, 0]]^m PDF[UniformDistribution[{0, d}], x], {x, 0, d}];

capDelta[a5_?NumericQ, a6_?NumericQ, d_?NumericQ] := expVal[a6, d, 1] expVal[a5, d, -1];

capBeta[a5_?NumericQ, a6_?NumericQ, d_?NumericQ] := expVal[a5, d, 1] expVal[a6, d, -1];

delta[a1_?NumericQ, a2_?NumericQ, a3_?NumericQ, a4_?NumericQ, d_?NumericQ] := 
         expVal[a4, d, 1] expVal[a2, d, -1] expVal[a3, d, -1] expVal[a1, d, 1];

beta[a1_?NumericQ, a2_?NumericQ, a3_?NumericQ, a4_?NumericQ, d_?NumericQ] := 
         expVal[a3, d, 2] expVal[a1, d, -2] expVal[a4, d, -1] expVal[a2, d, 1];

tbl[a1_?NumericQ,a2_?NumericQ,a3_?NumericQ, a4_?NumericQ,a5_?NumericQ, a6_?NumericQ] :=
   Transpose[{capDelta[a5, a6, #], capBeta[a5, a6, #], delta[a1, a2, a3, a4, #], 
              beta[a1, a2, a3, a4, #]} & /@  Table[i, {i, .1, 1, .1}]];

labels = {"Delta", "Beta", "delta", "beta"};

Using the definitions above:

Manipulate[lp = ListPlot[tbl[a1, a2, a3, a4, a5, a6], 
     Joined -> True, DataRange -> {.1, 1}, PerformanceGoal -> "Speed", ImageSize -> 400]; 
  linestyles = Cases[lp, {directive__, line_Line} :> {directive}, \[Infinity]]; 
  Row[{lp, Grid[Table[{Graphics[Append[linestyles[[i]], Line[{{-1, 0}, {1, 0}}]], 
        ImageSize -> 50, AspectRatio -> 1/10], labels[[i]]}, {i, 4}], 
        Spacings -> 2, Alignment -> Left]}], 
{{a1, .5, "a1"}, .1, 1, .1}, 
{{a2, .5, "a2"}, .1, 1, .1}, 
{{a3, .5, "a3"}, .1, 1, .1}, 
{{a4, .5, "a4"}, .1, 1, .1},
{{a5, .5, "a5"}, .1, 1, .1}, 
{{a6, .5, "a6"}, .1, 1, .1}]

screenshot:

enter image description here

share|improve this answer
    
Thanks again, this is great. On pen-and-paper though, I hoped to have Mathematica and esp. Manipulate as a nice environment to track how my observations would change with different errors, so it would be great if things would work numerically whatever I throw at it (non-uniform errors, correlation etc.). But thanks! –  László Jul 13 '12 at 20:20
add comment
up vote 0 down vote accepted

Manipulate of course works with CompoundExpression, sot the code below does what I was after. Thank you for all the help.

Manipulate[u={u1,u2,u3,u4,u5,u6};
  answers=MapThread[(Min[1,Max[#1+#2,0]])&,{u,{e1,e2,e3,e4,e5,e6}}];
  disc={answers[[5]],answers[[6]],answers[[3]]/answers[[1]],answers[[4]]/answers[[2]]};
  DB={disc[[2]]/disc[[1]],disc[[1]]^2/disc[[2]]};
  db={disc[[4]]/disc[[3]],disc[[3]]^2/disc[[4]]};
  deltasbetas[x_]:=NExpectation[{DB,db},{e1,e2,e3,e4,e5,e6}
   \[Distributed]UniformDistribution[Table[{0,x},{6}]]];
  Needs["PlotLegends`"];
  ListPlot[Transpose@Table[Flatten@deltasbetas[x],{x,0.01,0.51,0.05}],
   Joined->True,PlotLegend-> {"Delta","Beta","delta","beta"},
   LegendPosition->{1.1,-0.4},LegendShadow->None],
   {{u1,0.8},0,1,0.1},{{u2,0.7},0,1,0.1},{{u3,0.5},0,1,0.1},{{u4,0.4},0,1,0.1},
   {{u5,0.8},0,1,0.1},{{u6,0.7},0,1,0.1}]

The question is: Could it be any faster, like with some Compile or parallelization?

share|improve this answer
add comment

Still slow, but much faster than yours:

u = {0.8, 0.7, 0.5, 0.4, 0.8, 0.7};
answers = MapThread[(Min[1, Max[#1 + #2, 0]]) &, {u, {e1, e2, e3, e4, e5,  e6}}];
disc = {answers[[5]], answers[[6]], answers[[3]]/answers[[1]], answers[[4]]/answers[[2]]};
DB = {disc[[2]]/disc[[1]], disc[[1]]^2/disc[[2]]};
db = {disc[[4]]/disc[[3]], disc[[3]]^2/disc[[4]]};
deltasbetas[x_] := NExpectation[{disc, DB, db}, {e1, e2, e3, e4, e5, e6} \[Distributed] 
                   UniformDistribution[Table[{0, x}, {6}]]]

ListPlot[Transpose@Table[Flatten@deltasbetas[x], {x, .1, 1, .1}], Joined -> True]

enter image description here

share|improve this answer
    
Nice to keep it numeric, thanks! –  László Jul 13 '12 at 20:21
    
I got a $Failed (in version 8.0.4), as NIntegrate tried to e1 = x as a limit and could not. Does the code you pasted work for you, verbatim? Just curious, I suspect my mistake, as usual. Thanks! –  László Jul 13 '12 at 22:31
1  
@László It works here. Try starting a clean session –  belisarius Jul 13 '12 at 22:33
    
Strange, that indeed did the trick. Thanks! –  László Jul 13 '12 at 22:36
    
Now I only need the whole thing to respond to the numbers in u. :) (Well, and also labeling this many lines somehow.) –  László Jul 13 '12 at 23:16
show 1 more comment

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.