Advanced Search

janhardo

840 Reputation

12 Badges

11 years, 221 days
Contact janhardo

MaplePrimes Activity


These are replies submitted by janhardo

@Jean-Michel  The Position fun...

janhardo 840
January 20 2026
0 0

@Jean-Michel 

The Position function works even more extensively in Mathematica.

It depends on what you want to calculate.

Try to make this possible in Maple.

@Jean-Michel    Position := ...

janhardo 840
January 20 2026
0 0

@Jean-Michel 
 

Position := proc(expr, target, path := [])
    local i, res := []:

    if expr = target then
        return [path];
    elif type(expr, list) then
        for i to nops(expr) do
            res := [op(res), 
                    op(Position(expr[i], target, [op(path), i]))];
        od;
    elif type(expr, `+`) then
        # Voor sommen
        for i to nops(expr) do
            if type(op(i, expr), `*`) then
                # Voor producttermen in een som: match het hele product
                if has(op(i, expr), target) then
                    res := [op(res), [op(path), i]];
                end if;
            else
                # Voor andere termen: ga recursief
                res := [op(res),
                        op(Position(op(i, expr), target, [op(path), i]))];
            end if;
        od;
    elif type(expr, `*`) then
        # Voor producten buiten een som: ga recursief
        for i to nops(expr) do
            res := [op(res),
                    op(Position(op(i, expr), target, [op(path), i]))];
        od;
    fi;

    return res;
end proc:


# Test
M := [[1, a], a, [[3 + a]], 2*a - 1]:
Position(M, a);

matrices_of_matrices_mprimes_20-1-2026.mw

bvp)_bad_index_into_Matrix_mprimes_20-1-...

janhardo 840
January 20 2026
0 0

bvp)_bad_index_into_Matrix_mprimes_20-1-2026.mw

nm 11983  Isn't it fantast...

janhardo 840
January 20 2026
0 0

nm 11983 
Isn't it fantastic that Maple gives a dot for multiplication in expressions !

It clearly shows when there is multiplication and can be useful for some users, I think.

complete mixture pdf with modus...

janhardo 840
January 19 2026
0 0

@sand15 
"The mode  can be quite difficult to catch given the wide variety of pdf a mixture can have.
So this is a quite simple, likely not optimal, procedure to do this."

Can this procedure : FindGlobalMode_FP_Visual_RV_simple handle all mixtures ?

 
> 

restart;
with(plots):
with(Statistics):

FindGlobalMode_FP_Visual_RV_simple := proc(
    GRVs::list,           # Lijst van RandomVariable objecten
    pi_weights::list,     # Mengselgewichten
    N::posint := 10^4,
    tol::numeric := 1e-10,
    maxit::posint := 200
)
  local K, j;
  local starts, z, znew, i, k;
  local modes, vals, xmode;
  local f, Gmap, pi, mu_list, sigma_list;
  local sample_list, rnd_val, cum_prob, component_chosen;
  local xmin, xmax, max_std, overall_mean;
  local hist_plot, pdf_plot, mode_line, mode_point;
  local u1, u2, s, z_box;

  K := nops(GRVs);

  if nops(pi_weights) <> K then
    error "GRVs and pi_weights must have the same length";
  end if;

  # Gebruik Mean en StandardDeviation functies
  mu_list := [seq(Mean(GRVs[j]), j=1..K)];
  sigma_list := [seq(StandardDeviation(GRVs[j]), j=1..K)];  # sigma, niet sigma^2

  # Toon voor debugging
  print("Component means:", mu_list);
  print("Component standard deviations:", sigma_list);

  # Normalize weights
  pi := [seq(pi_weights[j]/add(pi_weights), j=1..K)];

  # --------------------------------------------------
  # 1. ANALYTISCHE PDF - via RandomVariable PDFs
  # --------------------------------------------------
  f := x -> add(pi[j]*PDF(GRVs[j], x), j=1..K);

  # --------------------------------------------------
  # 2. VASTE-PUNT MAP (gebruik sigma^2 = sigma_list[j]^2)
  # --------------------------------------------------
  Gmap := x ->
    add(pi[j]*mu_list[j]/(sigma_list[j]^2)
        *exp(-(x-mu_list[j])^2/(2*sigma_list[j]^2))/sqrt(2*Pi*sigma_list[j]^2),
        j=1..K)
    /
    add(pi[j]/(sigma_list[j]^2)
        *exp(-(x-mu_list[j])^2/(2*sigma_list[j]^2))/sqrt(2*Pi*sigma_list[j]^2),
        j=1..K);

  # --------------------------------------------------
  # 3. ZOEK GLOBALE MODUS
  # --------------------------------------------------
  starts := [op(mu_list), add(pi[j]*mu_list[j], j=1..K)];

  modes := [];
  vals  := [];

  for k to nops(starts) do
    z := evalf(starts[k]);

    for i to maxit do
      znew := evalf(Gmap(z));
      if abs(znew - z) < tol then
        break;
      end if;
      z := znew;
    end do;

    modes := [op(modes), z];
    vals  := [op(vals), f(z)];
  end do;

  xmode := modes[ ListTools:-Search(max(vals), vals) ];

  # --------------------------------------------------
  # 4. SAMPLING - DIRECTE BOX-MULLER METHODE
  # --------------------------------------------------
  sample_list := Vector(N, datatype=float[8]);
  
  # Initialiseer random seed
  randomize();
  
  for i from 1 to N do
    # Gebruik rand() voor uniform [0,1] om component te kiezen
    rnd_val := rand()/10^12;
    
    # Kies component
    cum_prob := 0;
    component_chosen := 1;
    for j from 1 to K do
      cum_prob := cum_prob + pi[j];
      if rnd_val <= cum_prob then
        component_chosen := j;
        break;
      end if;
    end do;
    
    # Genereer normaal verdeeld getal met Box-Muller
    s := 2;
    while s >= 1 do
      u1 := 2*rand()/10^12 - 1;
      u2 := 2*rand()/10^12 - 1;
      s := u1^2 + u2^2;
    end do;
    
    z_box := u1*sqrt(-2*ln(s)/s);  # Standaard normale variabele
    
    # Transformeer naar gekozen component
    sample_list[i] := mu_list[component_chosen] + sigma_list[component_chosen] * z_box;
  end do;

  sample_list := convert(sample_list, list);
  
  # --------------------------------------------------
  # 5. VISUALISATIE
  # --------------------------------------------------
  max_std := max(op(sigma_list));
  overall_mean := add(pi[j]*mu_list[j], j=1..K);
  
  xmin := min(overall_mean - 5*max_std, min(op(mu_list)) - 4*max_std);
  xmax := max(overall_mean + 5*max_std, max(op(mu_list)) + 4*max_std);

  hist_plot := Histogram(sample_list, color="LightGray", style=polygon);
  pdf_plot := plot(f(x), x=xmin..xmax, color=red, thickness=2);
  mode_line := plot([[xmode,0],[xmode,f(xmode)]], color=blue, linestyle=dash, thickness=2);
  mode_point := pointplot([[xmode, f(xmode)]], color=blue, symbol=solidcircle, symbolsize=15);

  print(display(hist_plot, pdf_plot, mode_line, mode_point,
    title=cat("Gaussian Mixture (K=", K, "): Global Mode = ", sprintf("%.6f", xmode)),
    labels=["x","density"],
    view=[xmin..xmax, DEFAULT]
  ));

  return xmode;

end proc:
> 

# Goed gescheiden componenten
G13 := RandomVariable(Normal(-4, 0.7)):
G14 := RandomVariable(Normal(0, 0.9)):
G15 := RandomVariable(Normal(4, 0.6)):

weights_sep := [0.33, 0.34, 0.33]:
mode_sep := FindGlobalMode_FP_Visual_RV_simple([G13, G14, G15], weights_sep, 5000);

"Component means:", [-4, 0, 4]

  

"Component standard deviations:", [.7, .9, .6]

  

 

3.999937260

 (1)
> 

# Goed gescheiden componenten
G13 := RandomVariable(Normal(-1, 0.7)):
G14 := RandomVariable(Normal(10, 0.9)):
G15 := RandomVariable(Normal(4, 0.6)):
# Goed gescheiden componenten
G16 := RandomVariable(Normal(-5, 0.7)):
G17 := RandomVariable(Normal(1, 0.9)):
G18 := RandomVariable(Normal(15, 0.6)):
# Goed gescheiden componenten
G19 := RandomVariable(Normal(-2, 0.7)):
G20 := RandomVariable(Normal(-1, 0.9)):
G21 := RandomVariable(Normal(13, 0.6)):
> 

weights_sep := [0.33, 0.34, 0.33, 0.35, 0.35, 0.30,0.4, 0.4, 0.4]:
mode_sep := FindGlobalMode_FP_Visual_RV_simple([G13, G14, G15, G16, G17, G18,G19, G20, G21], weights_sep, 5000);

"Component means:", [-1, 10, 4, -5, 1, 15, -2, -1, 13]

  

"Component standard deviations:", [.7, .9, .6, .7, .9, .6, .7, .9, .6]

  

 

-1.331489397

 (2)
 

NULL

Download mixture_mengsel_via_randon_variabele_invoer_normaal_verdeeld_pdf_mprimes_20-1-2026.mw

@sand15 Does it appear that a mixtu...

janhardo 840
January 19 2026
0 0

@sand15 ,thanks to make this clear.
Can it be that  my calculation is a special case ?

  > ....

janhardo 840
January 19 2026
0 0
> 

restart;
with(Statistics):

f := t -> w*PDF(Normal(0,a),t) + (1-w)*PDF(Normal(0,b),t);# mixture pdf

"maple.ini in user"

  

proc (t) options operator, arrow; w*Statistics:-PDF(Normal(0, a), t)+(1-w)*Statistics:-PDF(Normal(0, b), t) end proc

 (1)
> 

 

> 

assume(a > 0):
assume(b > 0):
assume(w >= 0, w <= 1):

f := t -> w*PDF(Normal(0,a),t) + (1-w)*PDF(Normal(0,b),t): # expextation

EZ := Int(t*f(t), t=-infinity..infinity)=simplify(int(t*f(t), t=-infinity..infinity));

Int(t*((1/2)*w*2^(1/2)*exp(-(1/2)*t^2/a^2)/(Pi^(1/2)*a)+(1/2)*(1-w)*2^(1/2)*exp(-(1/2)*t^2/b^2)/(Pi^(1/2)*b)), t = -infinity .. infinity) = 0

 (2)
> 

NULL
> 

EZ2 := simplify(int(t^2*f(t), t=-infinity..infinity));# second moment

(1-w)*b^2+a^2*w

 (3)
> 

NULL
> 

VarZ := simplify(EZ2 - EZ^2); # variance

(1/2)*(-(Int(-(a*(-1+w)*exp(-(1/2)*t^2/b^2)-w*exp(-(1/2)*t^2/a^2)*b)*t, t = -infinity .. infinity))^2+2*b^2*(-(-1+w)*b^2+a^2*w)*a^2*Pi)/(Pi*a^2*b^2) = (1-w)*b^2+a^2*w

 (4)
> 

NULL
> 

NULL

 

 
> 

 

> 

# Define mixture distribution with known mean and variance
mixtureDist := Distribution(
    PDF = unapply(w * PDF(Normal(0, a), t) + (1-w) * PDF(Normal(0, b), t), t),
    Mean = 0,
    Variance = w*a^2 + (1-w)*b^2,
    Conditions = [a > 0, b > 0, 0 <= w, w <= 1]
);

Z := RandomVariable(mixtureDist);

# Check properties
PDF(Z, t);                    # Shows mixture PDF
Mean(Z);                      # Returns 0
Variance(Z);                  # Returns w*a^2 + (1-w)*b^2

_m2773745101568

 

_R7

 

(1/2)*w*2^(1/2)*exp(-(1/2)*t^2/a^2)/(Pi^(1/2)*a)+(1/2)*(1-w)*2^(1/2)*exp(-(1/2)*t^2/b^2)/(Pi^(1/2)*b)

 

0

 

(1-w)*b^2+a^2*w

(1.1)
 

 

Download how_to_build_a_gausian_mixture_model_19-1-2026.mw

gausian_mixture_maple_primes_19-1-2026.mw

plotting_with_functions_mprimes_17-1-202...

janhardo 840
January 17 2026
0 0

plotting_with_functions_mprimes_18-1-2026.mw

correct Eq 2,9 ?All other derivations fr...

janhardo 840
January 15 2026
0 0




correct Eq 2,9 ?
All other derivations from Eq 2.9  in paper are wrong ?

 

table_for_ode_solution_25_cases_mprimes_...

janhardo 840
January 07 2026
0 0

table_for_ode_solution_25_cases_mprimes_7-1-2025.mw

Eq_14an_paper_afgeleid_mprimes_6-1-2026....

janhardo 840
January 06 2026
0 0

Eq_14an_paper_afgeleid_mprimes_6-1-2026.mw

Not using this : it are perodic function...

janhardo 840
January 01 2026
0 0

Not using this : it are periodic functions : t>=0,t<=2*Pi
They repeating themselves. 
 ,  ( )
generalized sin  function : JacobiSN(u, 0) = -1

sol := solve({cos(t) = 0, sin(t) = -1}, t, allsolutions);

works too for the generalized sin and cos function

sol := solve({
    JacobiCN(u, 0) = 0,
    JacobiSN(u, 0) = -1
}, u,allsolutions);

@C_R  Expr := 'int(1 / (sq...

janhardo 840
December 29 2025
0 0

@C_R

Expr := 'int(1 / (sqrt(-alpha*l^2 + 1) * sqrt(-alpha*l^2 * (x + 1)/2 + 1)), alpha = 0 .. z)';
EvaluatedExpr := value(Expr);
'int(1 / (sqrt(-alpha*l^2 + 1) * sqrt(-alpha*l^2 * (x + 1)/2 + 1)), alpha = 0 .. z)' = EvaluatedExpr;

@C_R This form is also acceptable ?...

janhardo 840
December 29 2025
0 0

@C_R 



This form is also acceptable ?, or the denominator must still  be 1 ?,although it is the same integral.
Expression II_subs is  chanced by a function.

What expression for the integrand is bet...

janhardo 840
December 29 2025
0 0

What expression for the integrand is better to handle ?

1 2 3 4 5 6 7 Last Page 2 of 82