sand15

787 Reputation

11 Badges

9 years, 185 days

MaplePrimes Activity


These are replies submitted by sand15

@Josolumoh 

In my previous reply, please change

 

MAX := max(rhs~(Tally(SampleB)));
B0  := eval(B(rhs~(NE)[]), x=0);
scale := MAX/B0;

 

by

 

 

MODE := sort(Tally(SampleB), key=(x->rhs(x)))[-1];  #return the mode
XMAX := lhs(MODE);
HMAX := rhs(MODE):

B0  := eval(B(rhs~(NE)[]), x=XMAX);
scale := HMAX/B0;


This will give a better scaling of the graph of the probability function if the mode is not at x=0.
(try for instance r=2, b=1, d=0.1)

 

Last but not least: if the value of Xmax is large enough for the probability that
x > Xmax is very small, then the scale is simply given by


scale := `number of samples`


PS: you can verify that the sampling method described in the previous reply gives as good results as Statistics:-Sampling in the case d=0 (negative binomial distribution).

Enjoy

@Josolumoh 


A spreviously quoted by Carl the case d=0 is obvious.
More precisely it corresponds to a negative binomial distribution with parameters r and b/(2+b)
(in Maple  NegativeBinomial(r, b/(2+b))  )

Carl also pointed the meaningless of the case d<0

I describe here the way to sample your  "pseudo distribution" B (which does not sum to 1 and then need to me normalized).
The last plot superimposes the probability function of B for r=3, b=2, d=1 (red line) and the histogram in discrete form.
I considered that only integers values of x have a sense.

For other valuers of r, b, d you will have to verify that the value of the variable Xmax is large enough (in this case the variable ShouldBe1 must have reached a "converged" value; to verify this increase the value of Xmax and look to what happened to ShouldBe1).

Let me know ih this suits you

 

restart

with(Statistics):

A := (x+r-1)!/(r-1)!/x! * p^r * (1-p)^x

factorial(x+r-1)*p^r*(1-p)^x/(factorial(r-1)*factorial(x))

(1)

# verify A is the probability function of a NegativeBinomial distribution of parameters (r, p)

ProbabilityFunction(NegativeBinomial(r, p), x) assuming x > 0:
convert(%, factorial);

factorial(x+r-1)*p^r*(1-p)^x/(factorial(r-1)*factorial(x))

(2)

A := unapply(A, p):

pi := (b/2) / (1+d*x+b/2);
ip := (1+d*x) / (1+d*x+b/2);

(1/2)*b/(1+d*x+(1/2)*b)

 

(d*x+1)/(1+d*x+(1/2)*b)

(3)

# verify ip = 1-pi

simplify(ip-(1-pi));

0

(4)

A(pi);

factorial(x+r-1)*((1/2)*b/(1+d*x+(1/2)*b))^r*(1-(1/2)*b/(1+d*x+(1/2)*b))^x/(factorial(r-1)*factorial(x))

(5)

B := unapply(A(pi) / (1+d*x), [r, b, d]);

proc (r, b, d) options operator, arrow; factorial(x+r-1)*((1/2)*b/(1+d*x+(1/2)*b))^r*(1-(1/2)*b/(1+d*x+(1/2)*b))^x/(factorial(r-1)*factorial(x)*(d*x+1)) end proc

(6)


Case d = 0

# When d=0 B(r, b, 0) is the probability function of a NegativeBinomial distribution of parameters (r, pi)
# where pi=b/(2+b)

simplify(subs(d=0, pi));

B(r, b, 0);
map(simplify, B(r, b, 0));

b/(2+b)

 

factorial(x+r-1)*((1/2)*b/(1+(1/2)*b))^r*(1-(1/2)*b/(1+(1/2)*b))^x/(factorial(r-1)*factorial(x))

 

GAMMA(x+r)*(b/(2+b))^r*2^x*(1/(2+b))^x/(GAMMA(r)*factorial(x))

(7)

# sampling B(x, b, 0) is then trivial:
#  1/ set the values of r and b, for instance 3 and 5
#  2/ use Sample(NegativeBinomial(3, 5), 10)


Case d > 0
Note that the case d < 0 has no  meaning (see lso Carl's answers)

# numerical example (NE)

NE := [r=3, b=2, d=1]:
B(rhs~(NE)[])

(1/2)*factorial(x+2)*(1-1/(x+2))^x/(factorial(x)*(x+2)^3*(x+1))

(8)

plot(B(rhs~(NE)[]), x=1..10)

 

# sequence of values of B(...) for integer values of x in [0, 100]

Xmax := 100:

PointwiseB := [ seq([x, B(rhs~(NE)[])], x in [$0..Xmax]) ]:

# given the plot above, assume the following quantity should be equal to 1
# if B was a true probability function

ShouldBe1 := add(op~(2, PointwiseB)):
evalf(%);

.2266874330

(9)

# Normalize B:

PointwiseB := [ seq([x, B(rhs~(NE)[]) / ShouldBe1], x in [$0..Xmax]) ]:
evalf(add(op~(2, PointwiseB)));

1.

(10)

# Draw the approximate cumulative function

c := CumulativeSum(op~(2, PointwiseB));

ApproximateCDF := zip((u, v) -> [u[1], v], PointwiseB, convert(c, list)):

plot(ApproximateCDF);    # not an exact representation... should be a stepwise function

c := Vector(4, {(1) = ` 1 .. 101 `*Array, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

 

# Sample ApproximateCDF:
#  1/ draw u from a Uniform(0, 1)
#  2/ find the

U := RandomVariable(Uniform(0, 1)):
ListTools:-BinaryPlace(op~(2, ApproximateCDF), Sample(U, 1)[1]);

1

(11)

SampleB := map(u -> ListTools:-BinaryPlace(op~(2, ApproximateCDF), u), Sample(U, 1000));

SampleB := Vector(4, {(1) = ` 1 .. 1000 `*Vector[row], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(12)

Tally(SampleB)

[0 = 550, 1 = 177, 2 = 92, 3 = 42, 4 = 24, 5 = 16, 6 = 13, 7 = 10, 9 = 5, 8 = 9, 11 = 7, 10 = 7, 13 = 3, 12 = 3, 15 = 3, 14 = 2, 19 = 1, 17 = 2, 22 = 4, 23 = 1, 20 = 1, 21 = 1, 27 = 1, 26 = 1, 31 = 3, 30 = 2, 36 = 1, 39 = 2, 33 = 3, 34 = 1, 44 = 2, 41 = 1, 40 = 2, 42 = 1, 63 = 1, 56 = 1, 73 = 1, 77 = 1, 67 = 1, 90 = 1, 94 = 1]

(13)

MAX := max(rhs~(Tally(SampleB)));
B0  := eval(B(rhs~(NE)[]), x=0);
scale := MAX/B0;

window := [0..20, default]:

plots:-display(
      Histogram(SampleB, frequencyscale=absolute, discrete, thickness=10),
      plot(scale*B(rhs~(NE)[]), x=0..1, color=red, thickness=3),
      plot(scale*B(rhs~(NE)[]), x=1..100, color=red, thickness=3),
      view=window
)

550

 

1/8

 

4400

 

 

plot(scale*B(rhs~(NE)[]), x=0..1, adaptive=true)

 

 

 


 

Download B_Distribution.mw

 

@tomleslie

 

Sorry, you're right, I forgot the "1/3" in the coefficient of x^3.
 

@tomleslie 

Not applicable: see tomleslie's answer

You got a wrong result with Optimization:-NLPSolve (?)

Given fexact is monotone decreasing over [0, 2] and its values are between [0, 1], given fapp is monotone increasing on [0, 2] with fapp(2) > 6789*2, it is obvious that the the maximum M of abs(fapp-fexact) is reached for x=2.

Then M = subs(x=2, fapp-fexact) =  - sin(2)/2 which is about  30543.54535 and not 2827.54535.

Using instead Optimization:-Maximize( abs(fapp-fexact), x=0..2) returns the good value.

More of this: your two calculations (L2 norm and L1 norm) give strangely the same result: could it be that Optimization:-Maximize(  abs(fapp-fexact), x=0..2, maximize) does not do what you were expecting for?

@tomleslie 

I apreciate these non-gurus explanations.
Thanks Tom

@acer 

Thanks acer for all these technical explanations.

For my understanding: I guess than you writting :-`+`(a, b) is a shortcut for Tolerances:-`+`(a, b) isn't it?
If so, is it not an unsafe way to proceed ifyou use two different packages containing different functions with the same name (I recently met this situation with packages simplex and geom which both contain functions named convexhull)?

 

 

i

@Carl Love

I had noticed the `+` returned by the command line  < with(Tolerances); > but I forgot that  `+` does not operate on a list but on a sequence; thus I forgot to use op or [] before running  `+`.
I like your (`+`@op) solution, easier to understand than using subsindets(Z, :-`+`, `+`);
Thanks

It seems to me that you wanted to answer the question EB1000 asked some hours ago?
Unfortunately, your answer turned out to be a question!

@tomleslie

Thanks Tom,
 I do not understand either why "evalr" seems necessary. In any case, thank you for keeping me from going in circles.

@acer 

Thank you acer.

As often there are several ways to tackle a problem and looking for the "optimal" solution seems to be an endless story.

By the way, I used myself the "subs({0="H", 1="T"},...)" method to return a solution in terms of "H" and "T.
Unfortunately writting
UseHardwareFloats := true:
V := floor~(SignalProcessing:-GenerateUniform(numThrows, 0, 2))
does not return integers making the {0="H", 1="T"} inoperative.
I tried the substitution {0.="H", 1.="T"}  but it didn't work neither.

I understand that floor, round, ceil do not return integers when the argument is a harware floating point number.
Is there a simple way to obtain, for instance, the integer part of harware floating point number?

@Carl Love 

 

Thanks Carl for your very documented answer.
I think I will  try to implement your suggestions, Not that the problem of calculating the probabilities of each sum deserves it, but rather out of personal satisfaction (in "practical" situation one would probably use Central Limit Theorems and the normal approximation when the number of dice exceeds a fex dozens).

By the way, I used in Q3 the package SignalProcessing and, from my student years I remember that convolution is generally be done, for effiency, by using FFT (2 directe + one inverse).
I do not know if SignalProcessing:-Convolution proceeds this way, but if it is the case, do you think its limitations you pointed out also impact FFT?

 

@radaar 

Have you thought about initializing the seed of the pseudo random generator?

This is a classical trap: executing twice the same worksheet generates exactly the same pseudo random numbers.To avoid this use:

Seed := randomize()

in a separate block justr after your top "restart" line

@Carl Love  Thanks for your remark.

It's rater difficult to compare the results produced by Q3 with the true ones. I did that for 3 and 4 dice and that was all good.
Obviously I was overconfident my assuming it would be same for higher values.

I suppose your observation definitely discredits Q3?

 

Maybe I'm wrong in my interpretation?
I understant A is a random variable with integer outcomes 1..6 and that each outcome has a different weight B

Would that be okay with you?
restart:
with(Statistics):
A := [$1..6];                   # A = [1, 2, 3, 4, 5, 6]
B := A *~ 10;                 # B = [10, 20, 30, 40, 50, 60]
P := B /~ add(B);           # P = [1/21, 2/21, ....6/21]
X := RandomVariable(ProbabilityTable(P));
N := 100;
Histogram(Sample(X, N), discrete);

# also:
plot(CDF(X, t), t=min(A)..max(A), discont=true)

Remark: if A was another set of outcomes, for instance A:=[1, 2, 3, 5, 8, 13], the plot of the histogram would be wrong (abscissas would be 1, 2, ..6). In this case you should correct the histogram plot by writting
Histogram(Sample(X, N), discrete, tickmarks=[[seq(i=A[i], i=1..6)], default])

Let me know if my answer is off the mark.

@tomleslie 

I agree: it's always better to try a simple approach.

Nevertheless this doesn't mean we can avoid some preliminary reflection.
When you write that the "penalty" function (I prefer personnaly the term of "objective" function) is  sum( (xmodel-xexp)^2+(ymodel-yexp)^2), you make implicit unwritten hypothesis.
This writting, while usual in least squares method, relies on underlyng assumptions:

  1. The disagreement between the model predictions (x*(t), y*(t)) and the experimental points (x(t), y(t)) come from random additive measurement errors.
  2. These errors are homoskedastic  https://en.wikipedia.org/wiki/Homoscedasticity
  3. The measurement errors on X and Y are iid (independent and identically distributed)
  4. There are gaussian
  5. There is no measurement error on t
  6. There is no measurement error on alpha neither on v__0
  7. ...I'm not even sure this list is exhaustive

Under these hypotheses, sum( (xmodel-xexp)^2+(ymodel-yexp)^2) represents, up to ad additive constant, the log-likelihood of the experimental data given
the predictions from the model.

If only one of the previous assumptions fails, this writting is false.

In the case (3) fails, it's easy to modify your penalty function in  sum( ((xmodel-xexp)/sX)^2+((ymodel-yexp)/sY)^2) where sX and sY are the standard deviations
of the measurement errors on X and Y.

If the homoscedasticity assumption fails the penalty function should be  sum( wX(tn)*(xmodel(tn)-xexp(tn))^2+wY(tn)*(ymodel(tn)-yexp(tn))^2) where  wX(tn) and
 wY(tn) are weights depending on t (the previous case being a particular situation where   wX(tn) = (1/sX)^2 and wY(tn) = (1/sY)^2 for each tn.

if (4) fails the writting is simply false: think of errors with uniform distributions!

if (5) fails the only correct way to account for measurement errors on t is a bayesian approach.
The often invoked "generalized least squares method" is not rigourous an would fail if (4) fails too.

I stop here.
Let me be clear, I'm not saying your work is meaningless. I just say that writtig the formula you used to express a "penalty" function has no sense if you
do not explicitely state the hypothesis that makes it reliable.
Too many people use statistical methods as if they were cooking recipes.
I thought this clarification was necessary in order Erik may form himself a correct idea about the true complexity of its problem, and thus about the most suited method to use in order to solve it

First 12 13 14 15 16 17 18 Last Page 14 of 24