## 6568 Reputation

8 years, 103 days

## Statistics:-LinearFit can do t...

Statistics:-LinearFit can do the job quite well. So there is no bug in Statistics:-LinearFit  (or at least not on this problem).

The issue comes from the fact that the problem is not relevant, per se, of the framework of Ordinary Least Squares (OLS) regression where we aim to assess the values of the parameters A and B of the model:

" Conditional expectation of Y given X=x = B*x+A "

The reason for this is that the slope is known.
In the attached file I show how a rewritting of the initial problem sets it back in the OLS framework. But this is neither immediate not natural.

However, the simplest method to assess the value of "a" comes from classical estimation results in OLS : the best unbiased estimator of "a" is just the mean of pts2 - 3*pts1

 > restart:
 > with(Statistics):
 > pts1 := Vector([3, 5, 21]);
 > pts2 := Vector([4, 6, 16])
 (1)

It couldn't be easier

 > # From standard results in Ordinary Least Squares (OLS) regression, # the Best Linear Unbiased Estimator (BLUE) of the intercept A # in the model: Expectation(Y | X) = B*X+A is just # BLUE(A) = Mean(Y) - BLUE(B)*Mean(X) # # Here BLUE(B) is supposed to be known and set to 3. # Then BLUE__A := Mean(pts2 - 3 *~ pts1); add (pts2 - 3 *~ pts1) / 3
 (2)

Using Statistics:-LinearFit

 > # The problem is not relevant of OLS but can be solved with LinearFit # (at the price of a few pirouettes) # # Rewrite the problem this way: #   1/ define Z as Y-3*X Z := (X, Y) -> Y - 3*X; #   2/ imagine you repeat 3 times the same experiment for a same value of the #      control variable "a" ("a" is the only regressor of the problem) #      Write Newpts2 the vector of the outcomes of these experiments Newpts2 := Z~(pts1, pts2); #   3/ define then the "design matrix" as an arbitrary constant vector, #      for instance < 1, 1, 1> Newpts1 := Vector([1, 1, 1]); #   4/ apply LinearFit as usual, remembering that the only regressor is "a" L := LinearFit([a], Newpts1, Newpts2, a); #   5/ Substitute "a" by the value it takes in the vector Newpts1 A := subs(a=1, L); convert(%, rational);
 (3)
 >
 >

## Maybe this?...

Simulation of a "true" Galton board is very difficult.
You have, in all rigor, to simulate the drop of a sphere (a disc in 2D) shocking many fixed cylindrical pins (discs in 2D).
The trajectory of a disc (assuming the problem is 2D) is described by a succession of free flights submitted to the gravitational field (rather simple to code), each of them ending when the falling disc hurts a fixed pin (another disc).
Determining the impact point is a difficult problem.
More of this you have to describe the physics of the shock (elastic or not), account for possible slipping and rolling...
Last but not least a true Galton board has a limited width an the pins are confined in a kind of funnel whose walls modify the fall at at boundaries of the board.

What I propose you is a simpler simulation based on a random walk on a "triangular lattice" (Pascal's lattice).
Starting from the position (i, j) a "ball" can move either to position (i+1, j+1) or (i+1, j-1) with equal probabilities. The process is repeated from the initial position (0, 0) a given number of times.
Think to this to some kind of "discrete Galton board".

Look to the animation in the attached file.
If it suits you I will give you some explanation about it... if required

WATCH OUT: I faced difficulties to load the file with the animations, so I commented the corresponding line (red written on yellow background).
To run the animation, once the plot after the command Fallouts(100, 10) has been displayed, just clik on it and run the animation by using the keys in the upper left of your Maple worksheet.

It's likely that this code can be written in clever and more performing way.

 > restart:

see for instance https://en.wikipedia.org/wiki/Random_walk

 > with(Statistics): with(plots): with(plottools):
 > randomize():
 > R := rand(0 .. 1): Fallouts := proc(NP, N)    local  Fallout:    local  GaltonBoard, depart, t, n, c, T, GBH;    global cs, GB:          # fallout of a single "ball"      Fallout := proc(np)       local GBH, depart, t, n, c, T;       global cs, GB:       GBH := copy(GB):       depart := [0, 0];       t      := depart:       for n from 1 to N do          depart := depart +~ [1, (2*R()-1)];          t      := t, depart       end do:          t  := [t, [depart[1]+1, depart[2]]];       c  := unapply(CurveFitting:-Spline(op~(1, t), op~(2, t), x, degree=1), x);       T  := transform((x,y) -> [y+N+1, x]):          if np > 1 then          GBH := display(GBH, T(Histogram([cs], discrete, thickness=5, frequencyscale=absolute)) )       end if;              cs := cs, depart[2];        plots[animate](plot, [c(x), x=0..A, color=red, thickness=2], A=0..N+1, frames=N+2, background=GBH);        end proc:    GaltonBoard := [seq(seq([i, j], j in [seq(-i..i, 2)]), i=0..N)]:    GB := plot(GaltonBoard, style=point, symbol=solidcircle, symbolsize=10, scaling=constrained, gridlines=true, color=blue):    cs := NULL:    plots[animate]( Fallout, [np], np=1..NP, frames=NP); end proc:
 > global cs: # drop 100 "balls" on a "discrete Galton Board" made oh 10 layers of pins # Fallouts(100, 10);   # TO BE UNCOMMENTED
 > cs
 (1)
 > # Histogram([cs])
 >

## Maybe you could try this   ...

Maybe you could try this

 > restart
 > with(Statistics):
 > XG := h -> RandomVariable(Normal(0, h))
 (1)
 > MaxIt := 1000
 (2)
 > randomize():
 > # Be careful: h must be strictly positive h := evalf([seq(1-k/MaxIt, k=0..MaxIt-1)]): h[1],h[-1];
 (3)
 > # do you really need to print the results ? # # for i from 1 to MaxIt do #    printf("%a\n%a\n", Sample(XG(h[i]), 5), i); # end do:
 > # Maybe this could suit you? # Stage 1 ; generate a random matrix by sampling Normal(0, 1) M := Sample(XG(1), [MaxIt, 5], method = [envelope, range = -1 .. 1]): # Watch out: stage 2 is computationnaly expensive. #            It is presented here to help understanding the method # # Stage 2: rescale M according to the standard deviations 1, 0.999, ...0.001 S := LinearAlgebra:-DiagonalMatrix(h): # Stage 3: do the product S.M Res := S.M
 (4)
 > # The same thing in a more efficient way M   := CodeTools:-Usage( Sample(XG(1), [MaxIt, 5], method = [envelope, range = -1 .. 1]) ): Res := CodeTools:-Usage( <  seq(M[n, ..] *~ h[n], n=1..MaxIt) > );
 memory used=0.73MiB, alloc change=0 bytes, cpu time=35.00ms, real time=11.00ms, gc time=0ns memory used=0.60MiB, alloc change=0 bytes, cpu time=9.00ms, real time=3.00ms, gc time=0ns
 (5)
 > for i from 1 to MaxIt do   # printf("%3d : %a\n", i, [entries(Res[i,..], nolist)]); end do:
 >

## tomleslie has provide a clever and fast ...

To complete their propositions, here is a time and memory comparison of different ways to answer yout problem

 > restart:
 > with(Statistics):   # the most natural (?) and one of the fastest if you want 0 and 1 outputs B := RandomVariable(Bernoulli(1/2)):
 > numThrows:= 10^5:
 > CodeTools:-Usage( Sample(B, numThrows) ):
 memory used=0.87MiB, alloc change=0 bytes, cpu time=7.00ms, real time=7.00ms, gc time=0ns
 > restart;
 > with(StringTools) : # a simple way to get "H" and "T" ouputs
 > numThrows:= 10^5: CodeTools:-Usage( Random(numThrows, "HT" ) ):                     # a string of the form "HHT..." CodeTools:-Usage( Explode(Random(numThrows, "HT" )) ):            # a list of the form [ "H", "H", "T", ... CodeTools:-Usage( Vector(Explode(Random(numThrows, "HT" ))) ):    # a vector (here fot comparison with tomleslie's solution below)
 memory used=154.93KiB, alloc change=0 bytes, cpu time=193.00ms, real time=193.00ms, gc time=0ns memory used=0.86MiB, alloc change=0 bytes, cpu time=192.00ms, real time=193.00ms, gc time=0ns memory used=22.27MiB, alloc change=0 bytes, cpu time=383.00ms, real time=384.00ms, gc time=41.38ms
 > restart;
 > # tomleslie  ... probably the fastest (?) to deliver "H" and "T" ouputs numThrows:= 10^5:
 > r:= rand(1..2):
 > CodeTools:-Usage(  Vector[row](numThrows, i->r()) ): CodeTools:-Usage(  Vector[row](numThrows, i->["H","T"][r()]) ):
 memory used=0.81MiB, alloc change=0 bytes, cpu time=33.00ms, real time=33.00ms, gc time=0ns memory used=0.77MiB, alloc change=0 bytes, cpu time=38.00ms, real time=38.00ms, gc time=0ns
 > restart;
 > with(LinearAlgebra): # advantage: there exist a very fast way to tansform 0 & 1 into "H" & "T"
 > numThrows:= 10^5:
 > V := CodeTools:-Usage(  RandomVector[row](numThrows,generator=rand(0..1)) ):  # outputs are 0 or 1 W := CodeTools:-Usage(  subs({0="H", 1="T"}, V) ):                            # outputs are "H" or "T"
 memory used=4.70MiB, alloc change=2.00MiB, cpu time=50.00ms, real time=51.00ms, gc time=5.63ms memory used=0.76MiB, alloc change=0 bytes, cpu time=3.00ms, real time=3.00ms, gc time=0ns
 > restart;
 > with(SignalProcessing):  # comparable to Statistics
 > numThrows:= 10^5:
 > V := CodeTools:-Usage( floor~(GenerateUniform(numThrows, 0, 2)) ):  # the fastest (?), use Hfloats (default) UseHardwareFloats := false: CodeTools:-Usage( floor~(GenerateUniform(numThrows, 0, 2)) ):  # don't use Hfloats
 memory used=1.79MiB, alloc change=0 bytes, cpu time=7.00ms, real time=8.00ms, gc time=0ns memory used=203.74MiB, alloc change=47.02MiB, cpu time=1.58s, real time=1.59s, gc time=88.87ms
 > restart:
 > with(combinat):  # just for fun
 > numThrows:= 10^5:
 > CodeTools:-Usage( randperm(["T"\$numThrows, "H"\$numThrows])[1..numThrows] ):
 memory used=20.67MiB, alloc change=4.58MiB, cpu time=161.00ms, real time=162.00ms, gc time=61.65ms
 >

## After some digressions unrelated to your...

After some digressions unrelated to your initial question, here is my definitive contribution

final.mw

 > restart:
 > with(Statistics):
 > Seed := randomize():
 > alias(RV=RandomVariable):
 > # Number of dice N := 3;
 (1)
 > # These two lines can be replaced by any algorithm amid those presented in this # same thread (just choose the one you want) S := add(RV(DiscreteUniform(1, 6)), k=1..N): P := Probability~(S =~ [\$N..6*N])
 (2)
 > # Define a discrete RV with outcomes of probabilities given by P # # Watch out: by default (I do not know if and how this could be changed), this RV has outcomes 1..5*N+1 # This raises an issue that will be fixed when the resuklts are ploted (the outcomes are N..6*N) Dice := RV(ProbabilityTable(P))
 (3)
 > # Sample "Dice" as any other RV M := 100: X := Sample(Dice, M):
 > h := Histogram(X, minbins=5*N+1): # The next operation is aimed to superimpose the exact mss function (aka P) and its # empirical estimation h. # # I transform h into a list that can be plot with ColumnGraph. # Watch out: #      Maple 2018: use op~(2, op~(1, [op(1..-2, op(1, h))] )); #      Maple 2015: use map(u -> print(u[-1][2]), [op(1..-2, op(1, h))]); interface(version); H := map(u -> u[-1][2], [op(1..-2, op(1, h))]): ColumnGraph(     [H, P],     legend=["Sampled", "Exact"],     color=[red, yellow],     background=gray,     title=cat("Sample size = ", M),     tickmarks=[[seq(i=i+N, i=0..5*N)], default],     size=[600, 400] )
 >

## @tomleslie  A few remarks: Use&...

@tomleslie

A few remarks:

1. Use Quantile(X, k+i+j, numeric) for a faster generation:
CodeTools:-Usage(Matrix(100, 10, (i,j) -> Quantile(X, (i-1)/100+(j-1)/1000, numeric))):
memory used=14.25MiB, alloc change=0 bytes, cpu time=158.00ms, real time=159.00ms, gc time=0ns
CodeTools:-Usage(Matrix(100, 10, (i,j) -> Quantile(X, (i-1)/100+(j-1)/1000))):
memory used=430.05MiB, alloc change=0 bytes, cpu time=6.28s, real time=5.90s, gc time=555.98ms

Note the CPU time (in this particular case) can be halved: Quantile(X, p) = -Quantile(X, 1-p)

2. Avoid computing Quantile(X, 0) which is obviously -infinity and not -15.681...

Here is a version which uses DocumentTools for a smarter rendering

QuantileTable.mw

## As you wrote it, your maplet can...

As you wrote it, your maplet can't work: Evaluate needs a function as argument, not a sequence of instructions.
Here is the correct way to code it in a simple case (no "for a in .... end do" loop ... but you will be able to generalize with from Carl's answer).

PS: the code I wrote conforms to your's, but it could be improved (remark first that there is no need to click on the "Display function" function button to visualize the MathML code in the MMLV field). It's also better to define "Input function  f(x)=" within a Label component, ...

 > restart:
 > with(Maplets[Elements]): with(Student[Calculus1]):
 > f := proc() local g, TL2: g := proc(func, slp)   local a;   a := Tangent(func,x=slp):   Maplets:-Tools:-Set('MMLV1'=a) end proc: TL2:=Maplet([   [     "Input function  f(x)=", TextField['func'](20)   ],   [     Button("Display function", Action(Evaluate('MMLV'='MathML[Export](func)')))   ],   MathMLViewer['MMLV'](),   [     "Input slope", TextField['slp'](10)   ],     [     Button("Tangent line", Action( Evaluate('function'=g, Argument(func), Argument(slp))))   ],      MathMLViewer['MMLV1']()    ]): Maplets[Display](TL2): end proc: f()
 >

You had forgotten that the projected dx*dy surface wasn't the surface itself:

 > restart;
 > with(VectorCalculus):
 > SetCoordinates(cartesian[x, y, z]):
 > R := sqrt(2): s := x^2+y^2+z^2-R^2; u := VectorField([x, 1, z]);
 (1)
 > n := N/sqrt(add(N[k]^2, k = 1 .. 3)); u . n; un := algsubs(x^2+y^2+z^2=R^2, %): un := simplify(un)
 (2)
 > intr1 := solve(subs(z = 0, s), y);
 > intr2 := solve(subs(z = 0, y = 0, s), x);
 (3)
 > # forgot that the projected dx*dy surface is not the surface irself l     := sqrt(x^2+y^2); h     := sqrt(2-x^2-y^2); phi   := arctan(h/l); scale := simplify(1/sin(phi)) assuming x^2+y^2 > 0; # surface of the sphere R1 := 2*int(int(scale, y = intr1[2] .. intr1[1]), x = intr2[2] .. intr2[1], numeric); evalf(4*Pi*R^2);
 (4)
 > R1 := 2*Student:-MultivariateCalculus:-MultiInt(scale*un, y = intr1[2] .. intr1[1], x = intr2[2] .. intr2[1])
 (5)
 > evalf(%);
 (6)
 > du := Divergence(u)
 (7)
 > intr1 := solve(s, z);
 > intr2 := solve(subs(z = 0, s), y);
 > intr3 := solve(subs(z = 0, y = 0, s));
 > R2 := int(int(int(du, z = intr1[2] .. intr1[1]), y = intr2[2] .. intr2[1]), x = intr3[2] .. intr3[1]);
 > evalf(%)
 (8)
 >

## In your particular case it's far mor...

In your particular case it's far more easier to verify the Gass theorem in spherical coordinates

I'll give a closer look to your code later

 > restart;
 > with(IntegrationTools): with(VectorCalculus): SetCoordinates(cartesian[x, y, z]):
 > MySphere := x^2+y^2+z^2=2; u := VectorField([x, 1, z]);
 (1)
 > n := N/sqrt(add(N[k]^2, k = 1 .. 3)); n := subs(MySphere, n);
 (2)
 > un := DotProduct(u, n); du := Divergence(u);
 (3)
 (4)
 > J := (r, phi) -> r^2*sin(phi); int(int(int(du*J(r, phi), theta=0..2*Pi), phi=0..Pi), r=0..SphereRadius); evalf(%);
 (5)
 (6)
 >

## Not the solution, just a way to derive t...

Not the solution, just a way to derive the equations by applying differential operators.
I think you have an ambiguity about the pressure (you declare p bar as a constant but it has a priori different values at z=-1 and z=+1 (unless pa is null)

For the solution: I'm not familiar with pdsolve and PDETools, but I'm afraid that you will have to do the things by hand.
If I'm not mistaken this is a classical poiseuille flow with uniform velocity at the inlet (parabolic profile expected at the outlet) and I remember solving it by hand long ago at school. Maybe you will have to reproduce explicitely the same modus operandi in Maple?

 > restart:
 > with(Student[VectorCalculus]):
 > SetCoordinates(cylindrical[r,theta, z]):
 > alias(u=u(r, 'theta', z)); alias(v=v(r, 'theta', z)); alias(p=p(r, 'theta', z));
 (1)
 > V := VectorField(< v, 0, u >); Divergence(V);
 (2)
 > Continuity_equation := eval(expand(%));
 (3)
 > Laplacian(V): Diffusion_term := beta^2 *~ map(expand, eval(%));
 (4)
 > Pressure_term := 1 *~ Gradient(p)
 (5)
 (6)
 > C := :-Vector[column](3, i -> Advection_term[i] = - Pressure_term[i] + Diffusion_term[i])[[1, 3]]
 (7)
 > NS := [Continuity_equation, entries(C, nolist)]
 (8)
 > # You can say "p bar is a constant and { p(z=-1)=0 & p(z=-1)=0 & p(z=-+1)=p__a } unless p__a=0 too   BC := [ :-D[1](u)(0, 'theta', z) = 0, v(0, 'theta', z) = 0, u(1, 'theta', z) = 0, v(1, 'theta', z) = p+alpha ]: seq(print(BC[i]), i=1..numelems(BC))
 (9)
 >

## Replace your line   Diam; by  ...

by                           Diam := simplify(Diam);

Nevertheless I would be tempted to say that Maple should simplify automatically.

Just rewrite it

PS:Personnaly I never use this document mode input that I find very capricious and subjected to hide characters mistakenly typed

## I'm not sure I understand you perfec...

I'm not sure I understand you perfectly.
What I think is: you have a discrete RV "X" with mass function fX(x)=1/(x+1) for x = 0, ..., k, where k is some given number (>0).

If It is that, a way to proceed is to use a ProbabilityTable Distribution

 > restart:
 > with(Statistics):
 > # Discrete distribution from a probability table # The sum of all masses must be equal to 1. P := k -> [seq(1/(x+1), x=0..k)] /~ add(1/(x+1), x=0..k) ; # example: K := 5: N := 10^3: P(K); X := RandomVariable(ProbabilityTable(P(K))): S := Sample(X, N): S[1..10]; Histogram(S);
 >

## Hi,  If your only concern is ...

Hi,

If your only concern is to compare the model and the experiment, and if this experiment is a "not too large" collection of triples (t, x, y) , then the option output=Array[ ...], where ... is the list of the experimental t values, is probably the simpler solution.

The answer is very simple: there is unfortunately no good way to proceed!
Computing "the Root Mean Square of all (xexperiment(t)-xmodel(t), yexperiment(t)-ymodel(t))" seems very reasonable.
But it can be a very poor solution for some problems, specifically if x or y are of different natures, vary on very different scales and so on.
Luckily it doesn't seem to be the case here.

Another point you must be aware of: if the (t, x, y) experimental values come from measurements, there are probably entailed by "errors". Suppose that x and y are perfectly measured, but that t is not.
For instance the "true" t is T=1, and the measured value is t=1.1: then comparing the experiment at time t=1.1 to the model at time t=1.1 doesn't make sense.
A (the?) correct approach in this case is to use a bayesian framework where one infer the true value of T given its measure t.

Maybe you could be interested in "fitting" your model on the experiments? For instance to assess the value of A ...?
You also mention that "The problem is however that even if the (x,y) curves are rather close, if may not be a good model".
Are you concerned with aapossible "model error"?
Here again, for these two questions, the correct answer relies to the "bayesian calibration of computer model" issue.

Let me know if you are interested by taht.

To conclude: try your  "Root Mean Square of all (xexperiment(t)-xmodel(t), yexperiment(t)-ymodel(t)", examine carefully the residues (the model-to-experiment distance at each time t) and if all looks good, then it's done

interface(rtablesize=max(Dimensions(MAT))):
MAT;

## Assuming the DEFINITION below suits you,...

Assuming the DEFINITION below suits you, you will find some elements in the attached file.

DEFINITION:

• Let E the ellipse defined by the equation (x/a)^2+(y/b)^2=1
• For each point P on E one define s(P) as the measure of the angle (counter clockwise counted) between Ox and OM
• Let A and B two points on E such that
• s(A) > =  0
• s(B) > s(A)
• s(B) <= 2*Pi
• Let M any point such that s(A) <= s(M) <= s(B)
• Let x(A) the abscissa of A and x(B) the one of B

Then the mean value of the distance d(O, M) from O=(0, 0) to M as M moves from A to B is defined by the integral of d(O, M) between x(A) and x(B), divided by x(B)-x(A).

In file Ellipse2.mw you will find a procedure 'md' which computes the mean value of d(O, M). It takes 4 arguments:

1. b
2. alpha = s(A)
3. beta = s(B)

'md' handles (if I'm not mistaken) the 3 cases:

1. alpha < Pi and beta < Pi
2. alpha < Pi and beta > Pi
3. alpha > Pi and beta > Pi

Before that you will find a few lines that treat more formally the first situation alpha < Pi and beta < Pi

As claimed in this file, I was not capable to obtain a general formal solution for any values of a and b (that's why I arbitrarily wrote a:=2 and b:=1 at the top at the file).
The search of a general formal solution requires using assumptions about a, b, p, q (see the commands in the file) and the 'int' procedure of Maple is not very comfortable with 'assume' or 'assuming' statements (someone here will probably correct me on this point)

Ellipse2.mw

 First 49 50 51 52 53 54 55 Page 51 of 56
﻿