MAINT: stats.multinomial: use dtype-dependent tolerance

[mdhaber]

Reference issue

gh-22585

What does this implement/fix?

gh-22585 adjusted the condition under which p is normalized in scipy.stats.multinomial. multinomial has converted p to float64 internally since it was introduced, but users have been allowed to pass in input with other dtypes, and the change did not take this into account. This compensates by making the tolerance dtype-dependent.

[mdhaber]

Suggested change

	eps = np.finfo(np.result_type(np.asarray(p), np.float32)).eps * 10
	eps = np.finfo(np.result_type(np.asarray(p), np.float16)).eps * 10

maybe?

[ev-br]

can this be just `eps = xp.finfo(xp_default_dtype(xp)).eps * 10 # or other wiggle factor`

This will not work currently since xp_default_dtype only account for torch defaulting to f32;
Assuming that the failure mode is that jax returns f32 array even when explicitly asked for an f64, xp_default_dtype should grow some logic to account for JAX idiosyncrasy.

EDIT: no, what I wrote above makes no sense: it's not related to JAX at all, it's a pure float32 thing #22585 (comment)

Then the right thing to do seems to be just eps = np.finfo(asarray(p)).dtype).eps * 10

[mdhaber]

There are tests that are passing integers like [0, 1] or something, so not without other changes.

[ev-br]

Ok then, that figures. Then indeed, there's no way around result_type. Sorry for the noise

[mdhaber]

@ev-br marked this as ready given #22585 (comment). I suppose we could also require that input be float64 if that conversion is going to be done anyway, but I'd rather not.

[ev-br]

LGTM, merged. Thanks Matt