Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 356 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

In Maple, it's a one-liner, requiring no local variables:

ProxySort:= (S::list,P::list)-> S[sort(P, 'output'= 'permutation')]:

The expression sqrt(x)/x is unstable in Maple. As soon as it's  created (whether by directly typing it or by rationalize), it gets automatically simplified to 1/sqrt(x)

Like this:

MyMatrix:= apply~(Matrix((2,3), symbol= elem), t);

I cannot get LinearFit to fit a single constant parameter. It seems that that case is explicitly forbidden. But it can be easily done with plain minimize:


f:= x-> 3*x+a:
S:= minimize(add((f~(pts1)-pts2)^~2), location);
plot(
   [<pts1|pts2>, eval(f(x), indets([S], identical(a)=anything))], 
   x= (min..max)(pts), 
   style= [point,line], symbolsize= 20
);

Which gives a = -61/3, same as Kitonum's method.

Use

simplify(evalc([(()->"Re")=Re, (()->"Im")=Im](Ec))) assuming positive;

Hint: Substitute epsilon= 1/upsilon, then simplify, then consider "by hand" the limit as upsilon approaches +infinity. I think that this is easier to see than the limit as epsilon approaches 0 from the right. Likewise, limit(..., epsilon= 0, left) is equivalent to limit(..., upsilon= -infinity).

And, of course, to make your problem more understandable to Maple, you must replace [ ] with ( ) and 0.5 with 1/2 (as VV said).

Here is an extendable template that should let you answer (in a natural way if you know statistical language) any such question for any discrete distribution with integer support. All results come directly from applying simplify (with assuming) and sum to the pmfs of the distribution. Each distribution is represented by a Maple container structure called a Record, which contains fields for the pmf, support, and assumptions on the parameters.

All parameters can be specified as any algebraic expressions (including constants); they need not be simple names. The `if` expressions in the Records adjust the assumptions appropriately.

The table (actually a DataFrame) at the end actually looks like a table when viewed in a worksheet. The MaplePrimes renderer squashes it here.
 

restart
:

Simplify:= (e,D)-> (simplify(e) assuming D:-assum[])
:

#Discrete expected value:
ExV:= proc(f, D::record(pmf, assum::set, support::range))
local x;
   Simplify(sum(f(x)*D:-pmf(x), x= D:-support), D);  
end proc
:

#
#Some abstract statistical parameters
#(feel free to add your own):
#

#Discrete mean:
Mean:= D-> Simplify(ExV(x-> x, D), D)
:

#Nth central moment:
Mom:= proc(n::nonnegint, D)
local x;
   Simplify(ExV(unapply((x-Mean(D))^n, x), D), D)
end proc
:

#Discrete variance:
Variance:= D-> Simplify(Mom(2,D), D)
:

#Sum-to-one check:
Unity:= D-> Simplify(ExV(1,D), D)
:

Skewness:= D-> Simplify(Mom(3,D)/Variance(D)^(3/2), D)
:

Kurtosis:= D-> Simplify(Mom(4,D)/Variance(D)^2, D)
:

CDF:= proc(D)
local X,n;
   unapply(Simplify(sum(D:-pmf(x), x= op(1,D:-support)..X), D), X)
end proc
:

#
#Some (parameterized) discrete distributions
#(feel free to add your own):
#

Bin:= (n,p)-> Record(
   'name'= 'Binomial'(n,p),
   'pmf'= (x-> binomial(n,x)*p^x*(1-p)^(n-x)),
   'assum'= {
      `if`(n::constant, [][], n::nonnegint),
      `if`(p::constant, [], [p >= 0, p <= 1])[]
   },
   'support'= 0..n
):

DiscrUni:= (a,b)-> Record(
   'name'= 'DiscreteUniform'(a,b),
   'pmf'= (x-> 1/(b-a+1)),
   'assum'= {
      `if`(a::constant, [][], a::integer),
      `if`(b::constant, [][], b::integer),
      `if`(andmap(type, [a,b], constant), [][],  a < b)
   },
   'support'= a..b
):

Poiss:= lambda-> Record(
   'name'= 'Poisson'(lambda),
   'pmf'= (x-> lambda^x/exp(lambda)/x!),
   'assum'= `if`(lambda::constant, {}, {lambda > 0}),
   'support'= 0..infinity
):

NegBin:= p-> Record(
   'name'= 'NegativeBinomial'(p),
   'pmf'= (x-> (1-p)^x*p),
   'assum'= `if`(p::constant, {}, {p >= 0, p <= 1}),
   'support'= 0..infinity
):

Table:= (DL::list(record), FL::list(procedure))-> DataFrame(
   Matrix((nops(DL),nops(FL)), (i,j)-> FL[j](DL[i])),
   columns= FL,
   rows= (D-> D:-name)~(DL)
):

Table(
   [Bin(n,p), DiscrUni(a,b), Poiss(lambda), NegBin(p)],
   [Unity, Mean, Variance, Skewness, Kurtosis]
);

DataFrame(Matrix(4, 5, {(1, 1) = 1, (1, 2) = p*n, (1, 3) = -p*n*(-1+p), (1, 4) = (1-2*p)/(sqrt(n)*sqrt(p)*sqrt(1-p)), (1, 5) = ((3*n-6)*p^2+(-3*n+6)*p-1)/(p*n*(-1+p)), (2, 1) = 1, (2, 2) = (1/2)*a+(1/2)*b, (2, 3) = (-(1/12)*b+(1/12)*a)*(-b+a-2), (2, 4) = 0, (2, 5) = (1/5)*(9*a^2+(-18*b-18)*a+9*b^2+18*b-12)/((-b+a)*(-b+a-2)), (3, 1) = 1, (3, 2) = lambda, (3, 3) = lambda, (3, 4) = 1/sqrt(lambda), (3, 5) = (3*lambda+1)/lambda, (4, 1) = 1, (4, 2) = (1-p)/p, (4, 3) = (1-p)/p^2, (4, 4) = -(p-2)/sqrt(1-p), (4, 5) = (-p^2+9*p-9)/(-1+p)}), rows = [Binomial(n, p), DiscreteUniform(a, b), Poisson(lambda), NegativeBinomial(p)], columns = [Unity, Mean, Variance, Skewness, Kurtosis])

 


 

Download DiscreteDist.mw

Joe has pointed out, essentially, that's it's often erroneous and usually suspicious to refer to a for loop's index variable outside the loop. The practice does have a few valid uses, such as determining whether the loop ended due to exceeding the maximum index or due to some other reason.

You have another major error: "Product of the eigenvalues" does not mean the dot product of the vector of eigenvalues with itself; rather, it means the (scalar) product of the three eigenvalues. You can do this with command mul, like this:

mul(N[i])

This is something like what your professor has in mind:

with(LinearAlgebra):
N:= 8:  n:= 3:  r:= true:
to N do
   M:= RandomMatrix(n,n);
   if Determinant(M) <> simplify(expand(mul(Eigenvalues(M)))) then 
      r:= false;
      break
   end if
end do: 
r;

I prefer this:

macro(LA=LinearAlgebra):
N:= 8:  n:= 3:  r:= true:
to N do
   if (radnormal@(LA:-Determinant - mul@LA:-Eigenvalues))(LA:-RandomMatrix((n,n))) <> 0 then
      r:= false;
      break
   fi
od:
r; 

I prefer it primarily because I trust radnormal when applied to a zero-equivalent radical expression more than I trust simplify@expand when applied to an integer-equivalent radical expression.

If is list (or set or sequence) of functions (or procedures) and is a sequence of arguments that would work for any of them, then simply F(S) is what you want. Note that this uses neither map nor an elementwise operator with ~. Example:

[iquo, irem](7,5)

or

(iquo, irem)(7,5)

I often use this as (min, max)(S). So, the number of arguments doesn't matter.  In particular, it could be one argument.

If and are lists of the same length, where are arguments for the first position and for the second, then do

F~(X, Y)

Unlike zip, this idea can be extended to any number of arguments. Any of the argument lists can be replaced by a single element, as long as that element is not itself a list, set, table, or rtable. The take-away lesson from all this is that there's no syntactic difference between a single function (or procedure) and a list, set, or sequence of them (an explicit sequence needs to be in parentheses).

Your code

map(eval~, [f(x), g(x)], x=~ [p,q,t])

can be (and should be) replaced by 

map(map, [f, g], [p, q, t])

because it's not necessarily possible to evaluate every function at symbolic x, and even if it's possible, it's not necessarily efficient.

The following procedure detects and reports cycles of exact algebraic numbers:

CycleDetectAlgnum:= proc(f, x0::radalgnum, Max::posint:= 999)
local x:= evala(Normal(x0)), R_inds:= table([x=0]), R_vals:= table([0=x]), k, j; 
   for k to Max do
      x:= evala(Normal(f(x)));
      if assigned(R_inds[x]) then return seq(R_vals[j], j= R_inds[x]..k-1) fi; 
      R_inds[x]:= k;  R_vals[k]:= x
   od;
   return
end proc
:   

It's usage is exactly like the other's:

CycleDetectAlgnum(y-> 4*y*(1-y), (5-sqrt(5))/8);

If you're only interested in detecting and reporting cycles rather than analyzing algebraically why they happen, then the following code---which does it all in hardware complex floats---will be much faster than either iterating the map or using a closed form from rsolve.

CycleDetectInner:= proc(
   C::Array(datatype= complex[8], order= C_order),
   k::integer[4],
   x::complex[8],
   eps::float[8]
)::integer[4];
option autocompile;
local j::integer[4];
   for j to k do
      if abs((C[j]-x)/`if`(abs(x) <= eps, 1, x)) <= eps then return j fi
   od;
   0
end proc
:
CycleDetect:= proc(
   f::procedure, 
   x0::complexcons, 
   Max::posint:= 999, 
   eps::{float,positive}:= evalhf(DBL_EPSILON)
)
option hfloat;
local x:= evalhf(x0), C:= Array(1..1, [x], datatype= complex[8], order= C_order), k, p; 
   for k to Max do
      x:= f(x);
      p:= CycleDetectInner(C, k, x, eps);
      if p > 0 then return C[p..] fi;
      C(k+1):= x
   od;
   return
end proc
:   

Example usage:

CycleDetect(y-> 4*y*(1-y), (5-sqrt(5))/8);

You may want to apply Re~ to the results to remove the zero imaginary parts. Any returned result is a 1-D Array containing a detected cycle. A return of NULL means that Max iterations were performed without a repeated value. 

The reason that add doesn't work as expected is that Tolerances doesn't overload it, whereas it does overload `+`. So, if you replace add with (`+`@op), you'll get your expected results.

Another option is post-processing of the partially evaluated results returned by add:

subsindets(Z, :-`+`, `+`);

The second `+` in that command is properly Tolerances:-`+`, but under the auspices of with(Tolerances), simply `+`  is sufficient.

I have two comments, unrelated to each other:

1. Regarding debug:

  1. The debug command cannot help you find syntax errors that prevent the procedure from even "compiling", which is the case here. I believe that the same is true for any language. To find such errors, pay attention to where your cursor is placed immediately after you get the error message. The error is often (but not always) within the same line.
  2. While you are typing code in Maple's GUI, and your cursor is immediately after a bracketing character (parenthesis, square bracket, etc.), a faint "shadow cursor" will appear on the matching bracket. Most errors of unbalanced brackets can be found this way.
  3. If you enter code in a Code Edit Region (available from the Insert menu), further syntax problems will be (subtly) highlighted as you type them.
     

2.  Regarding the time complexity of your algorithm:

The isolve command can find your solutions much faster (in general, as or -> infinity) than your brute-force double loop. Your procedure could be:

Fract:= proc(P::posint, Q::posint)  
local p,q;
   ((`/`@op)~)~(  #Convert integer pairs to fractions.
      select(  #Filter solutions to ranges 1..P-1 and 1..Q-1.
         type, 
         #Convert solution pairs to pairs of pairs:   
         map2(eval, [[p,q], [P-p, Q-q]], [isolve((P-p)*q-P*(Q-q) = 1)]),
         [[posint$2]$2]  #a pair of pairs of positive integers
      )
   )[]
end proc;  

 

Use sum to sum an infinite series, a finite indefinite series, or a finite definite series with a large number terms when it's suspected that  some symbolic tricks will work.

Your sum is 

sum(n/2^k, k= 0..infinity) + 1;

The resut is 2*n+1, not the 2*n-1 that you report. This is true for any n; it needn't be a power of 2.

Note that the 1 in your sum is a separate term; it doesn't fit the pattern of the other terms.

Not many infinite series can be handled by rsolve, but this one is especially simple since a closed form for the partial sums can be found. I'd use

RS:= rsolve({S(0)=1, T(1)= n, T(x)= T(x-1)/2, S(x)= S(x-1)+T(x)}, {T(x), S(x)});

where is the recurrence for the terms and S is the recurrence for the partial sums. To get the sum of the infinite series, take the limit of S(x) as approaches infinity:

limit(eval(S(x), RS), x= infinity);

Maple only supports 64-bit hardware floats[*1]. Software floats can be virtually any length, limited only by availlable memory and the exponent being a 1-word integer.

The following procedure converts between hardware floats and 16-character hex strings. If the argument can be represented as a hardware float but hasn't been, then a conversion to hardware float is performed. When a string is returned as a float, it is returned in a one-element Array. This is to prevent unintentional conversion to software float.

`hf<->hex`:= proc(x::{And(string, 16 &under length), Or(hfloat, realcons)}) 
   option hfloat; 
   `if`(x::string, hfarray(sscanf(x, "%y")), sprintf("%Y", evalhf(x)))
end proc:

#Example of use
Hx:= `hf<->hex`(sin(1));
                    Hx := "3FEAED548F090CEE"
Hf:= `hf<->hex`(Hx);
                   Hf := [0.841470984807897]

A useful excercise would be an additional procedure to extract the mantissa and exponent from the hex string (which'll require subdivision of at least one hexit, perhaps more). If you need, I'll do it. For software floats, this extraction is a trivial application of command op.

[*1]The term double precision is used colloquially throughout the English-speaking technical world for 64-bit hardware floats, and they are formally described by the IEEE-754 standard. It'd probably be best for you to either avoid or update your technical usage of the adjectives singledouble, etc. See Wikipedia IEEE 754.

 

First 151 152 153 154 155 156 157 Last Page 153 of 395