text:
Cremation is the hottest trend in the funeral industry
Jan. 22, 2013, 7:10 PM UTC
By Tyler Mathisen, Anchor
OK, forgive me. I know that headline was sort of cheesy. But when you've spent the past few months working on a documentary on the business of death, you go for the cheap laugh whenever you can.
I've been to casket-manufacturing companies, funeral homes, mortuary science schools and the annual convention of the National Funeral Directors Association. You don't know fun until you've partied with 6,000 funeral professionals.
They were a sliver of the more than 130,000 people who work in what's come to be known as the "death care industry" – a $17- billion-a-year business in America. And what's the fastest growing part of it? Cremation.
Every year in America, 2.5 million people die. In 2011, the last year for which numbers are available, 42 percent were cremated, according to the funeral directors association. That's double the rate of just 15 years ago. In some states, largely in the West, the cremation rate tops 70 percent. In Washington, it's 72 percent; in Nevada, almost 74 percent. (The lowest rate of cremation, in case you need a great pick-up line, is Mississippi's, at 15.7 percent.)
Why Rates are Rising: Religion, Family and Cost
So why the big jump in cremations? There are lots of reasons. One is the softening of the Catholic church's views of the practice. For centuries – until 1963, in fact – the church outlawed it. The church's laws still express a preference for burial. But the outright ban is a thing of the past and now, under some circumstances, bishops can permit a funeral mass with cremated remains present.
Another reason for the rise in cremations is the decline in nuclear families. As more Americans live far from hometowns and parents, and as family burial plots have waned in popularity and accessibility, millions have turned to cremation as a practical and cost-effective way to care for a loved one's remains.
But the main reason, as you might expect, is cost. Cremation is cheaper than burial. The average cost of a funeral today is about $6,500, including the typical $2,000-or-more cost of a casket. Add a burial vault, and the average jumps to around $7,700. A cremation, by contrast, typically costs a third of those amounts, or less. In a tough economy like the current one, cost counts – a lot.
Money Matters, Even in DeathWhat you pay depends in part on where you live and which additional products or services you buy from the funeral home or crematory. Prices tend to be higher in densely populated urban areas. And you will pay more – in some cases almost as much as a full burial with casket – if you contract with a funeral home for such ancillary services as hearses, visitations and viewings or memorial services prior to cremation.
In that case, you will be charged a few hundred dollars for – yes, I'm serious – a rental casket with removable liner. In most states, whether you have a viewing or not, you will have to buy a cremation container, usually wood fiber or cardboard, for $100 or so. It is burned with the body.
If you keep it simple, though, the average cost of a cremation, including a basic memorial service, runs about $1,600. If you go for a so-called direct cremation, without a memorial service, the cost can drop well below $1,000. Go online and you can find prices as low as $600 or so.
In the End: What to Do With the Ashes
Of course, the total cost doesn't stop with the cremation itself. You've got to do something with the 4 to 6 pounds of mostly carbon ash that remain after the body has been incinerated for several hours at 1,600 to 1,800 degrees Fahrenheit. Typically that means buying an urn. Cost: from under $100 to well over $1,000. Was your beloved Uncle Pete a golfer? You can get The Masters Golf Cremation urn online for $218. It won't improve your putting, or his.
Then there's the cost of interring the cremation urn. You don't have to buy a burial space for the remains. You can keep them in an urn on your mantel. Bad idea if you have cats or ball-throwing small children. Or you can scatter the ashes or have them mixed with concrete and dropped into the ocean to form a memorial reef. One vendor in Florida, Eternal Reefs, charges $3,000 to $7,000 for that. Says owner George Frankel, without a trace of irony: "I think a sea change is what we're seeing" in the funeral business.
If you decide to place your loved one's urn and ashes in a burial space, vault or columbarium, figure on spending $1,000 or so -- maybe more, depending on how fancy the cemetery or memorial garden is. You may pay extra, too, if you choose to have a niche with a view. But who's counting? Or looking, for that matter?
"Death: It's a Living" premieres Thursday Jan. 31 at 9 p.m. EST/PST on CNBC.
Tyler Mathisen, Anchor


text:
A weird week in underage drinking?
This weekly roundup in underage drinking news leaves us with more questions than answers.
Justin Bieber plays beer pong?
Maybe our little boy singer is growing up and making mistakes, but is everything reported about him the truth? Last week, we heard he had too much pasta and milk before performing and ended up vomiting onstage. Any parent could have warned him about that. Now, photos from last year – unearthed by TMZ yesterday – purportedly show him playing beer pong. Really? While we have zero verifiable information if he drinks too much beer or milk, the awkward looking photograph of a 17-year-old Bieber, wearing a hat, a half-zip fleece and a backpack only shows him looking at a bunch of red plastic cups on a ping-pong table while a half-naked unidentified man in the foreground holds the single beer shown. While it appears a beer pong game was in progress, nothing in the evidence demonstrates Bieber was playing.
Is a parent who hosts a 15-year-old's birthday party responsible when one party crasher stabs another?
In November 2006, a 15-year-old's birthday party with a half-dozen invited guests turned into a free-for-all when it was crashed by up to fifty uninvited teens, including warring gang members. Although the mom, Elizabeth Cumba, intervened to stop a fight that began in her basement, even personally marching the young man off her property to send him home, her attempts to control the situation were fruitless. Hiram Colon Jr., 17, stepped back onto the property to continue the argument, was chased by six teens, who police identified as East Hartford gang members, onto a neighbor's property and was fatally stabbed.
According to Lawyers USA, Cumba had successfully defended herself against a wrongful death lawsuit by Colon's estate. Now, "a state appeals court just ordered a new trial, finding the jury in her first trial decided the case under the wrong premises liability standard." What that means is the Connecticut appeals court decided the original judge's instructions to the jury were overly narrow in defining the law.
What makes matters worse is that none of the six youths who assaulted Colon received jail time. According to judicial abuse blogger Steven G. Erickson, who followed the case closely in his Stark Raving Viking blog, they were all offered a plea deal. Meanwhile, Colon's mother is left to mourn a child who died too young and Cumba, the mom who tried to get Colon to leave before the situation disintegrated, heads back to court. Like many parental nightmares, this one is set on repeat.
Murphy, Pat. "Benchmarks: Social host faces liability for fatal stabbing." Lawyers USA 26 Sept. 2012. General OneFile. Web. 10 Oct. 2012.
Why does an Assistant US Attorney gets a second chance we wouldn't be offered?
Video Link: Channel 9WSYR – ABC
In Casenovia, NY, it seems that the social host law is so confusing, even Assistant US Attorney Craig Benedict didn't understand how to protect his family and himself from the 100 to 200 teens who took over the family's 25 acres at a graduation party the Benedicts threw for their son. Photos from the party available on Instagram and social media messages intercepted afterwards by the police, indicated the teens planned to drink and "party" fully. That was in July. Now, a special prosecutor assigned to the case said that the case against Craig Benedict "will likely be dismissed if he stays out of trouble." Special treatment or professional courtesy?
This entry was posted on October 8, 2012 by Helene M. Epstein in Alcohol Policy, Minimal Proof required, Parental Responsibility, Parental Rights, Social Host, Social Host Ordinances, Underage Drinking and tagged Beer pong, Connecticut, Justin Bieber, Legal drinking age, Local Laws, Local ordinance, New York, Parenting, Santa Clara County California, Social Host Law, Social Host Laws, Underage Drinking.
https://wp.me/pCa0F-Hc

redpajama_set_name: RedPajamaCommonCrawl


Sitting on a windy precipice on the coast of Ireland, Chelsea, Fred and I discuss our experience with a short term missions trip to Ireland!
Are you a Human Being or a Human Do-ing?
I've heard it before—Life is more about the journey than the destination. You've heard it, too? When I hear that statement, I nod my head in agreement, like I know. Uh huh.
Last week was one of those weeks. I was running from one thing to the next, pencil poised over the to-do list, anxious to check things off. A wedding loomed on Saturday, written on my calendar for months, staring at me every day that previous week. Don't get me wrong. We are good friends with the parents of the groom. I love hanging out with these people. I just wasn't in a good place. I was weary and preoccupied and simply worn out.
We went to the wedding, and it was like medicine to the soul. New love. Old friends. Crazy pictures. ☺ As I sat there watching everyone talking, loving, communicating around the table at the reception, it felt really good.
Life is a journey. Birth, marriage, baptism, missions trips, Bible studies, picnics, church services, concerts, plays, movies, vacations, dinners together and even the Super Bowl. Let's break it down even further: laughing together, crying together, praying together, sharing together, learning together, playing together, helping each other. We want to share these episodes together. It's what makes up life. Not the destination.
I know this to be true. But I get so caught up with crossing items off my to-do list. There's a quote that says "Am I more focused on the 'to do' than I am on what I'm called 'to be'? Since I read that, I've made it a point to "be" more. To help more. But I wrestle with my wanting to get things done. It's a constant struggle.
How often am I like that dad? I remembered that story almost daily when Kailey was growing up. Now that she's an adult and moved out of the house, it still applies. That's why I want to take more time out to relax and enjoy family and friends, pray, laugh or help others. I know very well how hard it is when life requires us to "do." I'm just thinking it would be cool if we could work it in our crazy schedule every once in awhile. The journey is what makes life wonderful. Doing things with each other, for each other, enjoying each other.


text:
Q: How do I debug a Rust executable that uses a shared library? My current crate (X) I wish to debug uses another library crate of mine (Y), which loads a C shared library with Rusts ffi. The executable from crate X work fine, I can even debug it fine using GDB. However, when I want to debug the executable in visual studio code with CodeLLDB I am hit with "/<dir>/target/debug/<exe name>: error while loading shared libraries: .so: cannot open shared object file: No such file or directory".
I am puzzled as the library exists, Rust can build the executable, and GDB can debug perfectly fine without errors, but CodeLLDB cannot. So, my question is, how can I debug this Rust project with CodeLLDB, as this provides a nicer debugging experience (GDB has no interpretation of Rust in anyway)?

A: I was sloppy it seems. I think I only added the library directory to my .bashrc.
I now added a symbolic link to /usr/lib and debugging in VS Code works as expected.

redpajama_set_name: RedPajamaStackExchange


text:
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Linq.Expressions;
using System.Reflection;

namespace BufferedDataReaderDotNet.Extensions
{
    internal static class TypeExtensions
    {
        private const BindingFlags DefaultBindingFlags = BindingFlags.Instance | BindingFlags.Public;

        private static readonly IDictionary<Type, Func<BinaryReader, object>> ReadFuncs =
            new Dictionary<Type, Func<BinaryReader, object>>
            {
                {typeof(byte[]), br => br.ReadBytes()},
                {typeof(DateTime), br => br.ReadDateTime()},
                {typeof(Guid), br => br.ReadGuid()},
                {typeof(TimeSpan), br => br.ReadTimeSpan()}
            };

        private static readonly IDictionary<Type, Action<object, BinaryWriter>> WriteActions =
            new Dictionary<Type, Action<object, BinaryWriter>>
            {
                {typeof(byte[]), (o, bw) => bw.WriteBytes((byte[]) o)},
                {typeof(DateTime), (o, bw) => bw.WriteDateTime((DateTime) o)},
                {typeof(Guid), (o, bw) => bw.WriteGuid((Guid) o)},
                {typeof(TimeSpan), (o, bw) => bw.WriteTimeSpan((TimeSpan) o)}
            };

        static TypeExtensions()
        {
            AddReadFuncs();
            AddWriteActions();
        }

        private static void AddReadFuncs()
        {
            var methodInfos = typeof(BinaryReader)
                .GetMethods(DefaultBindingFlags)
                .Where(mi => mi.Name.Length > 4 && mi.Name.StartsWith("Read"))
                .Where(mi => mi.GetParameters().Length == 0)
                .Where(mi => mi.ReturnType != typeof(void))
                .GroupBy(mi => mi.ReturnType)
                .Where(g => g.Count() == 1)
                .Select(g => g.First());

            foreach (var methodInfo in methodInfos)
            {
                var returnType = methodInfo.ReturnType;

                if (ReadFuncs.ContainsKey(returnType))
                    continue;

                var binaryReaderParameter = Expression.Parameter(typeof(BinaryReader), "binaryReader");
                var callExpression = Expression.Call(binaryReaderParameter, methodInfo);
                var convertExpression = Expression.Convert(callExpression, typeof(object));

                var lamdbaExpression = Expression.Lambda<Func<BinaryReader, object>>(
                    convertExpression, binaryReaderParameter);

                ReadFuncs.Add(returnType, lamdbaExpression.Compile());
            }
        }

        private static void AddWriteActions()
        {
            var methodInfos = typeof(BinaryWriter)
                .GetMethods(DefaultBindingFlags)
                .Where(mi => mi.Name == nameof(BinaryWriter.Write))
                .Select(mi => new {Parameters = mi.GetParameters(), MethodInfo = mi})
                .Where(x => x.Parameters.Length == 1)
                .ToDictionary(x => x.Parameters[0].ParameterType, x => x.MethodInfo);

            foreach (var keyValuePair in methodInfos)
            {
                var parameterType = keyValuePair.Key;
                var methodInfo = keyValuePair.Value;

                if (WriteActions.ContainsKey(parameterType))
                    continue;

                var valueParameter = Expression.Parameter(typeof(object), "value");
                var convertExpression = Expression.Convert(valueParameter, parameterType);
                var binaryWriterParameter = Expression.Parameter(typeof(BinaryWriter), "binaryWriter");
                var callExpression = Expression.Call(binaryWriterParameter, methodInfo, convertExpression);

                var lambdaExpression = Expression.Lambda<Action<object, BinaryWriter>>(
                    callExpression, valueParameter, binaryWriterParameter);

                WriteActions.Add(parameterType, lambdaExpression.Compile());
            }
        }

        public static Func<BinaryReader, object> GetReadFunc(
            this Type fieldType, BufferedDataOptions bufferedDataOptions)
        {
            Func<BinaryReader, object> readFunc;

            if (!bufferedDataOptions.ReadFuncs.TryGetValue(fieldType, out readFunc))
            {
                if (!ReadFuncs.TryGetValue(fieldType, out readFunc))
                {
                    // TODO: methodInfo is a misnomer
                    var methodInfo = fieldType.GetMethods(BindingFlags.Public | BindingFlags.Static)
                        .Where(mi => mi.Name == "Parse")
                        .Select(mi => new {Parameters = mi.GetParameters(), MethodInfo = mi})
                        .Where(x => x.Parameters.Length == 1)
                        // If there's any ambiguity we'll prefer Parse(String).
                        .OrderByDescending(x => x.Parameters[0].ParameterType == typeof(string))
                        .First(); // InvalidOperationException if we can't read.

                    var binaryReaderParameter = Expression.Parameter(typeof(BinaryReader), "binaryReader");
                    var readStringMethodInfo = typeof(BinaryReader).GetMethod(nameof(BinaryReader.ReadString));
                    var callReadExpression = Expression.Call(binaryReaderParameter, readStringMethodInfo);
                    // We're trying to get to methods like SqlGeography.Parse(SqlString)
                    // whose SqlString parameter has an implicit conversion from String.
                    //var parseParameter = Expression.Convert(callReadExpression, methodInfo.Parameters[0].ParameterType);
                    //var callParseExpression = Expression.Call(null, methodInfo.MethodInfo, parseParameter);
                    //var convertExpression = Expression.Convert(callParseExpression, typeof(object));

                    // Just testing...
                    var convertExpression = Expression.Convert(callReadExpression, typeof(object));

                    var lambdaExpression = Expression.Lambda<Func<BinaryReader, object>>(
                        convertExpression, binaryReaderParameter);

                    return lambdaExpression.Compile();
                }
            }

            return readFunc;
        }

        public static Action<object, BinaryWriter> GetWriteAction(
            this Type fieldType, BufferedDataOptions bufferedDataOptions)
        {
            Action<object, BinaryWriter> writeAction;

            if (!bufferedDataOptions.WriteActions.TryGetValue(fieldType, out writeAction))
            {
                if (!WriteActions.TryGetValue(fieldType, out writeAction))
                    writeAction = (o, bw) => bw.Write(o.ToString());
            }

            return writeAction;
        }
    }
}

redpajama_set_name: RedPajamaGithub


text:
\section{Introduction} \label{s:intro}

Let $\Omega^*$ denote the set of $m\times n$ binary matrices with row sums $\boldsymbol r=(r_1,\dotsc,r_m)$ and column sums $\boldsymbol c=(c_1,\dotsc,c_n)$, and let $\boldsymbol w=(w_{ij})\in[0,\infty)^{m\times n}$ be a given nonnegative matrix.  Define the distribution $P^*$ on $\{0,1\}^{m\times n}$ via
\begin{equation} P^*(\boldsymbol z) = \frac{1}{\kappa} \prod_{ij} w_{ij}^{z_{ij}} {\mathds{1}}\{\boldsymbol z\in \Omega^*\} , \quad \quad \kappa = \sum_{\boldsymbol z\in\Omega^*} \prod_{ij} w_{ij}^{z_{ij}} , \label{e} \end{equation}
where ${\mathds{1}}$ is the indicator function and where we assume $\kappa > 0$.  $P^*$ is the conditional distribution of an $m\times n$ array of independent Bernoulli random variables, say $\boldsymbol B=(B_{ij})$, with $\Prob(B_{ij}=1)=w_{ij}/(1+w_{ij})$ given the margins $\boldsymbol r$ and $\boldsymbol c$, where $\Prob$ denotes probability.
This paper describes an importance sampling algorithm that can be used for Monte Carlo approximation of probabilities and expectations under $P^*$ and also for Monte Carlo approximation of $\kappa$.  After a preprocessing step, sampling from our proposal distribution requires $O(md)$ operations per matrix, where $d=\sum_j c_j=\sum_i r_i$ is the total number of ones in the matrix.

We are not aware of any existing importance sampling algorithms that permit practical inference under $P^*$, although many special cases have been studied in the literature.  For example, if $\boldsymbol w\equiv 1$ then $P^*$ is the uniform distribution over zero-one tables with specified margins, or equivalently, the uniform distribution over bipartite graphs with specified degree sequence.  For square matrices, if $\boldsymbol w$ is identically one except with a zero diagonal, then $P^*$ corresponds to the uniform distribution over directed graphs with specified degrees.  And if $\boldsymbol r = \boldsymbol c\equiv 1$, then $P^*$ is a distribution over weighted permutation matrices and $\kappa$ is the permanent of $\boldsymbol w$.  Empirically, our algorithm outperforms all existing importance sampling algorithms in these special cases.  Although our algorithm works well for most examples arising in practice, performance depends on $\boldsymbol r$, $\boldsymbol c$, and $\boldsymbol w$.  Highly irregular margins or highly variable $\boldsymbol w$, particularly many zero entries in $\boldsymbol w$, tend to cause poor performance.

$P^*$ factors in such a way that we need only focus on the distribution, say $P$, of the first column.  The columns are sampled sequentially, with each successive column viewed as the first column of a smaller matrix with updated margins based on the previously sampled columns.  We decompose the structure of $P$ into margin constraints, combinatorial factors, and non-uniform weighting terms, combine approximations of these terms in a principled way, and then use a dynamic programming algorithm to exactly and efficiently sample from the resulting proposal distribution $Q$ for the first column.  Sequentially sampling columns in this way defines a proposal distribution $Q^*$ for the whole matrix.  This strategy for algorithm design works well for many similar problems, including symmetric matrices and nonnegative integer-valued matrices, each of which will be described elsewhere owing to space constraints.  It seems likely that the design principles used for our approach are applicable much more broadly.

\section{Motivating applications}

\subsection{Conditional inference for graphs and tables} \label{s:conditional}

Let $\boldsymbol B\in\{0,1\}^{m\times n}$ be a matrix of independent Bernoulli random variables with
\begin{equation} \textstyle \logit \Prob(B_{ij}=1) = \alpha_i + \beta_j + \sum_k \theta_k \xi_{kij} , \label{e:logit} \end{equation}
where $\boldsymbol\alpha$, $\boldsymbol\beta$, and $\boldsymbol\theta$ are parameters, perhaps with constraints to ensure identifiability, and $\boldsymbol \xi$ is a collection of observed covariates.  Models of this form arise, for example, in educational testing, where $B_{ij}$ indicates whether or not subject $i$ responded correctly to question $j$.  If $\boldsymbol\theta\equiv 0$, then the model reduces to the classical Rasch model \citep{rasch1960probabilistic,rasch1961general}.  Otherwise, it is an extension of the Rasch model to include item-specific covariate effects, such as each subject's prior exposure to the content being tested in each question.  This model is also a simple version of models used for the analysis of network data \citep[c.f.,][]{holland1981exponential,fienberg1985statistical,goldenberg2010survey} where $\boldsymbol B$ is the adjacency matrix of a directed graph, $\boldsymbol \alpha$ and $\boldsymbol\beta$ allow for degree heterogeneity, and $\boldsymbol \xi$ is a collection of edge-specific covariates.  For example, for social network data we might have that $B_{ij}$ indicates whether subject $i$ reported subject $j$ as a friend, $\alpha_i$ controls the relative propensity for subject $i$ to report friends, $\beta_j$ controls the relative propensity for subject $j$ to be reported as a friend, and $\xi_{kij}$ indicates whether the relationship between subject $i$ and $j$ is of type $k$.

In both of these examples, if $\boldsymbol\theta$ is the only parameter of interest, then the nuisance parameters $\boldsymbol\alpha$ and $\boldsymbol\beta$ complicate inference and can be removed by conditioning on the row and column sums of $\boldsymbol B$ \citep[e.g.,][]{cox1958regression,holland1981exponential,mehta1995exact,harrison2012conservative}.  Conditioning also results in inferential procedures that are robust to modeling assumptions, implicit in \eqref{e:logit}, about the distribution of the margins.  The resulting conditional model that drives inference is exactly $P^*=P^*_\theta$ with $w_{ij}=\exp(\sum_k \theta_k \xi_{kij})$, perhaps with the additional constraint that $w_{ii}=0$ in the case of network data.  The conditional model is a natural exponential family in $\boldsymbol \theta$ with no nuisance parameters, but with an intractable normalization constant $\kappa=\kappa_\theta$.  \citet[Example 4.2]{harrison2012conservative} provides details about an importance sampling approach to exact conditional inference for this model.  The example there is based on a preliminary version of the algorithm presented here.

Other approaches to conditional inference in this setting include exhaustive enumeration, such as the algorithms for conditional logistic regression in Stata \citep{stata} and LogXact \citep{logxact}, Markov chain Monte Carlo approaches, such as the {\tt elrm} R package \citep{zamar2007elrm}, and analytic approximations, such as the {\tt cond} R package \citep{cond,brazzale2008accurate}, none of which are practical for larger matrices and/or multivariate $\boldsymbol\theta$.  Approximation of $\kappa$ was considered in \citet{barvinok2010number}.

\subsection{The uniform distribution and model validation}

If $\boldsymbol w\equiv 1$, or more generally, if $w_{ij}=\exp(\alpha_i+\beta_j)$ for real-valued $\boldsymbol \alpha$ and $\boldsymbol \beta$, then $P^*$ is the uniform distribution on $\Omega^*$.  The uniform distribution can be used for testing if $\boldsymbol\theta\equiv 0$ in model \eqref{e:logit}, or equivalently, for model validation of \eqref{e:logit} specialized to the case of $\boldsymbol\theta\equiv 0$.  Prominent examples include goodness-of-fit tests for the Rasch model \citep[e.g.,][]{rasch1960probabilistic,rasch1961general,ponocny2001nonparametric,chen2005exact,Chen:Sequential:2005} and for random bipartite graphs and directed graphs without reciprocity \citep[e.g.,][]{wasserman1977random,holland1981exponential,snijders1991enumeration}.

The uniform distribution over $\Omega^*$ also plays a central role in testing for the presence of interactions in co-occurrence tables, particularly co-occurrence tables arising in ecology, where $B_{ij}$ indicates the existence of species $i$ in location $j$ \citep[e.g.,][]{connor1979assembly,snijders1991enumeration,gotelli2000null,Chen:Sequential:2005}.  In these contexts, the uniform distribution subject to the margin totals is taken as a null hypothesis of no interaction among species (without necessarily assuming model \eqref{e:logit}).  Our original motivation for developing these algorithms came from a similar problem in neuroscience where $B_{ij}$ indicated whether neuron $i$ produced an action potential in time bin $j$, and the uniform distribution was taken as a null hypothesis of a lack of interaction among the neurons.  This can be viewed as an example of conditional testing for multivariate binary time series.  In that example, $mn\approx 10^8$ and existing algorithms in the literature were not practical.

Besides the many statistical applications, when $\boldsymbol w\equiv 1$, the normalization constant $\kappa$ is the number of binary matrices with specified margins, a topic of enduring interest in theoretical computer science \citep[e.g.,][]{kannan1999simple,jerrum2004polynomial,bezakova2007sbc} and combinatorial approximation \citep[e.g.,][]{bekessy1972asymptotic,mckay1984amp,greenhill2006aes,canfield2008aed,barvinok2010matrices}.  Importance sampling algorithms for $P^*$ can be used to provide efficient approximations of $\kappa$ \citep{blanchet2006isa}.

Monte Carlo sampling algorithms for the uniform distribution have been developed by many authors \citep[e.g.,][]{besag1989gmc, mckay1990uniform,snijders1991enumeration,rao1996markov,Chen:Sequential:2005, blanchet2006isa, bezakova2007sbc, chen2007cit, verhelst2008ema, Bayati:Sequential:2009}.  The approach here was inspired by the importance sampling algorithm in \citet{Chen:Sequential:2005}, but provides a more principled method for algorithm design that leads to substantial improvements in the uniform case and that also extends to the non-uniform case.

\subsection{Permanents and permanental processes}

If $\boldsymbol w$ is square and $\boldsymbol r=\boldsymbol c\equiv 1$, that is, if $\Omega^*$ is the set of permutation matrices, then $\kappa$ is the permanent of $\boldsymbol w$, also of enduring interest in theoretical computer science \citep[e.g.,][]{valiant1979complexity,jerrum2004polynomial}.  A variety of generalizations of permanents and determinants can be expressed as $\kappa\mu$, where $\mu=\Exp(h(\boldsymbol Z))$ for some function $h$, where $\Exp$ denotes expected value, and where $\boldsymbol Z$ has distribution $P^*$ \citep[e.g.,][]{littlewood1950theory,vere1988generalization,vere1997alpha,diaconis2000immanants}.  In principle, the algorithms here could be used to approximate the value of any of these objects, but the practicality of this approach depends heavily on $h$.  For example, the $\alpha$-permanent \citep{vere1988generalization,vere1997alpha} is
\begin{equation} \text{per}_\alpha(\boldsymbol w) = \kappa \Exp(\alpha^{\text{cyc}(\boldsymbol Z)}) , \label{e:perm} \end{equation}
where $\alpha\in\mathbb{R}$ and $\text{cyc}(\boldsymbol z)$ is the number of disjoint cycles in the permutation corresponding to $\boldsymbol z$.  The case $\alpha=1$ corresponds to the permanent of $\boldsymbol w$, and the case $\alpha=-1$ corresponds to $(-1)^n$ times the determinant of $\boldsymbol w$.  Permanents and $\alpha$-permanents arise in probability, statistics, and statistical physics in connection to permanental processes and random fields \citep[e.g.,][]{macchi1975coincidence,diaconis2000immanants,shirai2003random,mccullagh2006permanental,kou2009approximating} and the distribution of order statistics \citep{vaughan1972permanent,bapat1989order}.  Our approach is often effective for approximating \eqref{e:perm} when $\alpha>0$ and $|\log\alpha|$ is small.

\section{Algorithm design}

\subsection{The target distribution for the first column} \label{s:P}

For a matrix $\boldsymbol z=(z_{ij})$ we use $\boldsymbol z^j$ to denote the $j$th column of $\boldsymbol z$, we use $\boldsymbol z^{j:k}$ to denote the submatrix formed from columns $j,\dotsc,k$, we use $\boldsymbol R(\boldsymbol z)=(R_i(\boldsymbol z))$ to denote the column vector of row sums defined by $R_i=\sum_j z_{ij}$, and we use $\boldsymbol C(\boldsymbol z)=(C_j(\boldsymbol z))$ to denote the row vector of column sums defined by $C_j=\sum_i z_{ij}$.  Fix the size of the matrix, $m\times n$, the weights $\boldsymbol w$, and the margins, $\boldsymbol r$ and $\boldsymbol c$, and let $\boldsymbol Z$ have distribution $P^*$ defined in \eqref{e} with $\Omega^*=\{\boldsymbol z\in\{0,1\}^{m\times n}:\boldsymbol R(\boldsymbol z)=\boldsymbol r,\ \boldsymbol C(\boldsymbol z)=\boldsymbol c\}$.

To sample from $P^*$ we need only design a generic algorithm (generic in $m,n,\boldsymbol r,\boldsymbol c,\boldsymbol w$) for sampling from the distribution $P$ of the first column, namely,
\[ P(\boldsymbol x)=\Prob(\boldsymbol Z^1=\boldsymbol x) . \]  The reason is that the conditional distribution of $\boldsymbol Z^{2:n}$ given $\boldsymbol Z^{1}$ has the same form as $P^*$ in \eqref{e}, but with different parameters.  The size of the matrix is now $m\times(n-1)$, the row sums are updated to $\boldsymbol r-\boldsymbol R(\boldsymbol Z^1)$, the column sums are updated to $\boldsymbol c^{2:n}$, and the weight matrix is updated to $\boldsymbol w^{2:n}$.  Once we have sampled $\boldsymbol Z^1$, then we can update these parameters and effectively start over, treating the second column of $\boldsymbol Z$ like it was the first column of the new, updated problem, and then continuing sequentially until we have sampled the entire matrix.  This is the same sequential strategy suggested by \citet{Chen:Sequential:2005}.  The supplementary material (located at the end of this document) contains a more detailed description of this column-wise factorization.

Henceforth, our target distribution is $P$, the distribution of $\boldsymbol Z^1$.  Let $\boldsymbol Y$ be a random matrix chosen uniformly over $\Omega^*$, let $\Omega$ denote the support of $\boldsymbol Y^1$, namely,
\[ \Omega = \{\boldsymbol x\in\{0,1\}^{m\times 1}:\boldsymbol x = \boldsymbol z^1, \boldsymbol z\in\Omega^*\} , \]
and for $\boldsymbol x\in\Omega$ define
\[ U(\boldsymbol x) = \Prob(\boldsymbol Y^1=\boldsymbol x) , \quad \quad V(\boldsymbol x)=\Exp\Bigl(\prod_{ij} w_{ij}^{Y_{ij}} \Bigl| \boldsymbol Y^1=\boldsymbol x\Bigr) . \]
 It is straightforward to verify that
\begin{equation}  P(\boldsymbol x) \propto U(\boldsymbol x)V(\boldsymbol x){\mathds{1}}\{\boldsymbol x\in\Omega\}
\label{e:Pparts}
\end{equation}
for $\boldsymbol x\in\{0,1\}^{m\times 1}$.
This factorization conceptually isolates the hard margin constraints, $\Omega$, the combinatorics, $U$, and the non-uniform weighting, $V$.  Although the separation is clearly artificial, it is useful to treat each of these factors separately when developing a proposal distribution.

\subsection{The proposal distribution for the first column}

Motivated by the factorization in \eqref{e:Pparts}, we consider proposal distributions for the first column of the form
\[ Q(\boldsymbol x) \propto \tilde U(\boldsymbol x) \tilde V(\boldsymbol x) {\mathds{1}}\{\boldsymbol x\in\tilde\Omega\} , \]
where $\tilde U$ and $\tilde V$ are approximations of $U$ and $V$, respectively, that factor according to
\[ \tilde U(\boldsymbol x) \propto \prod_{i=1}^m u_i^{x_i} , \quad \quad \tilde V(\boldsymbol x) \propto \prod_{i=1}^m v_i^{x_i} \]
for some $\boldsymbol u,\boldsymbol v\in[0,\infty)^{m\times 1}$, and where $\tilde\Omega$ is of the form
\begin{equation} \label{e:Omega-tilde} \tilde\Omega = \bigl\{ \boldsymbol x \in\{0,1\}^{m\times 1} : x_{\pi_i}\in\mathcal{A}_i, \  \textstyle \sum_{\ell=1}^i x_{\pi_\ell}\in\mathcal{B}_i , \ i=1,\dotsc,m \bigr\} \end{equation}
for some permutation $\boldsymbol \pi=(\pi_1,\dotsc,\pi_m)$ of $(1,\dotsc,m)$ and some subsets $\boldsymbol{\mathcal{A}}=\mathcal{A}_1\times\dotsb\times\mathcal{A}_m\subseteq\{0,1\}^m$ and $\boldsymbol{\mathcal{B}}=\mathcal{B}_1\times\dotsb\times\mathcal{B}_m\subseteq\{0,1,\dotsc,c_1\}^m$.  Combining these approximations creates a proposal distribution of the form
\begin{equation} \label{e:Q}
 Q(\boldsymbol x) \propto \prod_{i=1}^m u_i^{x_i}v_i^{x_i}{\mathds{1}}\bigl\{x_{\pi_i}\in \mathcal{A}_i, \ \textstyle \sum_{\ell=1}^i x_{\pi_\ell} \in \mathcal{B}_i \bigr\} \quad\quad\quad (\boldsymbol x\in\{0,1\}^{m\times 1}) .  \end{equation}
Any proposal distribution of this form permits fast, exact sampling and evaluation using $O(mc_1)$ operations; see Section \ref{s:DP}.  The challenge is to find easily computable choices of $\boldsymbol u$, $\boldsymbol v$, $\boldsymbol \pi$, $\boldsymbol{\mathcal{A}}$, and $\boldsymbol{\mathcal{B}}$ such that $Q$ is a good approximation to the target $P$.  Fortunately, this seems to be possible in many cases; see Section \ref{s:components}.

For importance sampling to work, the support of $Q$, which is a subset of $\tilde\Omega$, must contain the support of $P$, which is a subset of $\Omega$.  When $\boldsymbol w$ has no zero entries, we require $\boldsymbol u$ and $\boldsymbol v$ to be positive and we engineer $\tilde\Omega$ to exactly coincide with $\Omega$ so that $P$ and $Q$ both have support $\Omega$; see Section \ref{s:margin}.  When $\boldsymbol w$ does have zero entries, we modify $\boldsymbol v$ and $\tilde\Omega$ to exclude certain elements in $\Omega$, but only elements that are not in the support of $P$.  This ensures that the support of $Q$ contains the support of $P$, but the supports may no longer be identical.
In this case, if the importance sampling algorithm generates a column that is not in the support of $P$, then as it sequentially generates additional columns it will eventually try to create a $Q$ that is identically zero, indicating that no assignment of the current column simultaneously satisfies the margin constraints and has positive weight.  At this point the algorithm can assign an importance weight of zero and terminate.  Certain patterns of zero weights make the algorithm highly inefficient because the algorithm rarely terminates with a nonzero importance weight.  In the supplementary material we discuss alternative choices of $\boldsymbol v$ and $\tilde\Omega$ that are more efficient for certain patterns of particular interest, including the special case of zeros only on the diagonal.

\subsection{Efficient sampling and evaluation of the proposal} \label{s:DP}

Let $\boldsymbol X\in\{0,1\}^{m\times 1}$ have a distribution $Q$ that factors according to \eqref{e:Q} above for some $\boldsymbol u$, $\boldsymbol v$, $\boldsymbol\pi$, $\mathcal{A}$, and $\mathcal{B}$.  Define the permuted partial sums $\boldsymbol S\in\{0,\dotsc,c_1\}^{m\times 1}$ according to $S_i=\sum_{\ell=1}^i X_{\pi_\ell}$ for each $i$, and note that $\boldsymbol X$ and $\boldsymbol S$ are in bijective correspondence.  The distribution of $\boldsymbol S$ factors according to
\begin{equation} \Prob(\boldsymbol S=\boldsymbol s) \propto \prod_{i=1}^m h_i(s_{i-1},s_i) \label{e:S} \end{equation}
for $h_i(s_{i-1},s_i) = u_{\pi_i}^{s_i-s_{i-1}}v_{\pi_i}^{s_i-s_{i-1}}{\mathds{1}}\{s_i-s_{i-1}\in \mathcal{A}_i, \  s_i \in \mathcal{B}_i \}$,
where here and below we define $S_0= s_0= 0$ for notational convenience.

The factorization in \eqref{e:S} implies that $\boldsymbol S$ is a Markov chain.  If we were given the standard Markov chain representation
\begin{equation}  \Prob(\boldsymbol S=\boldsymbol s) = \prod_{i=1}^m \Prob(S_i=s_i\bigl|S_{i-1}=s_{i-1}) \label{e:Qpi} ,\end{equation}
then generating a random observation of $\boldsymbol S$ would be trivial.  It is known that dynamic programming can be used to convert from Gibbs random field representations like \eqref{e:S} into Bayesian network representations like \eqref{e:Qpi}; see, e.g., \citet{Frey:Graphical:1998}.  The next theorem, which is straightforward to verify \citep[cf.][]{harrison2009rate}, summarizes dynamic programming in this context.

\begin{theorem} \label{t:DP} Let $(S_0,S_1,\dotsc,S_m)$ be a sequence of random variables where each $S_i$ takes values in the finite set $D_i$ and where $D_0=\{0\}$.  Suppose there exists a sequence of functions $h_i:D_{i-1}\times D_i\mapsto[0,\infty)$ for $i=1,\dotsc,m$ such that the distribution of $(S_1,\dotsc,S_m)$ can be expressed as
\[  \Prob\bigl(S_1=s_1,\dotsc,S_m=s_m\bigr) \propto \prod_{i=1}^m h_i(s_{i-1},s_i) . \]   Recursively define $g_m(s_{m-1},s_m) = h_m(s_{m-1},s_m)$ and
\[ g_i(s_{i-1},s_i) = h_i(s_{i-1},s_i)\sum_{s_{i+1}\in D_{i+1}} g_{i+1}(s_i,s_{i+1}) \quad (i=1,\dotsc,m-1), \]
where each $g_i$ is defined over $D_{i-1}\times D_i$.  Then $S_0,\dotsc,S_m$ is a Markov chain and
\[ \Prob\bigl(S_i=s_i\bigl|S_{i-1}=s_{i-1}\bigr) = \frac{g_i(s_{i-1},s_i)}{\sum_{t\in D_{i}} g_i(s_{i-1},t)} \quad (i=1,\dotsc,m) . \]
\end{theorem}

In the present context, $S_i\in\mathcal{B}_i\subseteq\{0,\dotsc,c_1\}$ for each $i$, so the algorithm described in Theorem \ref{t:DP} for converting from \eqref{e:S} to \eqref{e:Qpi} requires at most $O(mc_1^2)$ operations.  In fact, since in the present situation we have $h_i(s_{i-1},s_i)=0$ for $s_i-s_{i-1}\not\in\{0,1\}$, implying the same for $g_i$, this yields an algorithm that requires $O(mc_1)$ operations.  Instead of representing all $(c_1+1)^2$ combinations of $(s_{i-1},s_i)$, we represent only the $2c_1+1$ feasible combinations.  Once the representation in \eqref{e:Qpi} is computed, generating a random observation $\boldsymbol X$ from $Q$ or evaluating $Q(\boldsymbol x)$ at any $\boldsymbol x$ takes $O(m)$ operations.

\section{Specification of components} \label{s:components}

\subsection{Margin constraints} \label{s:margin}

Here we discuss the construction of $\tilde\Omega$. In particular, the next theorem shows how to ensure that $\tilde\Omega=\Omega$ for easily computable choices of $\boldsymbol\pi$, $\boldsymbol{\mathcal{A}}$, and $\boldsymbol{\mathcal{B}}$.
\begin{theorem} \citep{Chen:Sequential:2005} \label{t:binary} Assume $\Omega^*\neq\emptyset$. Choose $\boldsymbol\pi$ so that $r_{\pi_1}\geq \dotsb\geq r_{\pi_m}$.  For each $i=1,\dotsc,m$, define
\[ \mathcal{A}_i = \begin{cases}  \{0\} & (r_{\pi_i}=0);
 \\ \{0,1\} & (0 < r_{\pi_i} < n); \\
 \{1\} & (r_{\pi_i}=n), \end{cases} \quad \quad \quad
\mathcal{B}_i = \begin{cases} \{\max\{0,b_i\},\dotsc,c_1\} & (i< m) ; \\
\{ c_1\} & (i = m) , \end{cases}
\]
for $b_i = \sum_{\ell=1}^i (r_{\pi_\ell} -  \sum_{j=2}^{n}{\mathds{1}}\{c_j \geq \ell\})$.
Define $\tilde\Omega$ according to \eqref{e:Omega-tilde}.  Then $\tilde\Omega=\Omega$.
\end{theorem}

It is instructive to see how these choices of $\boldsymbol\pi$, $\boldsymbol{\mathcal{A}}$, and $\boldsymbol{\mathcal{B}}$ ensure that $\tilde\Omega\supseteq\Omega$, which is the primary requirement for importance sampling.  The Gale--Ryser conditions \citep{gale1957tfn,ryser1957cpm} state that there is a binary matrix with margins $\boldsymbol r\in\{0,\dotsc,n\}^{m\times 1}$ and $\boldsymbol c\in\{0,\dotsc,m\}^{1\times n}$ if and only if $\sum_i r_i=\sum_j c_j$ and
\[  \sum_{\ell=1}^i r_{\pi_\ell} \leq \sum_{\ell=1}^i \sum_j {\mathds{1}}\{c_j\geq \ell\} \quad \text{for all $i=1,\dotsc,m-1$}, \]
where the permutation $\boldsymbol \pi$ is chosen so that $r_{\pi_1}\geq\dotsb\geq r_{\pi_m}$.
It is straightforward to see that $\boldsymbol x \in\{0,1\}^{m\times 1}$ will be in $\Omega$ exactly when there is a way to fill out the remaining $n-1$ columns of the binary matrix that obey the updated margins after accounting for $\boldsymbol x$.  In other words, $\boldsymbol x\in\Omega$ exactly when $\boldsymbol r-\boldsymbol x$ and $\boldsymbol c^{2:n}$ satisfy the Gale--Ryser conditions for the margins of an $m\times(n-1)$ binary matrix.

The set $\boldsymbol{\mathcal{A}}$ is chosen so that $\boldsymbol x\in\boldsymbol{\mathcal{A}}$ if and only if $\boldsymbol r-\boldsymbol x\in\{0,\dotsc,n-1\}^{m\times 1}$.  The set  $\mathcal{B}_m$ is chosen so that $\sum_{i=1}^m x_i\in\mathcal{B}_m$ if and only if $\sum_{i=1}^m  r_i-x_i = \sum_{j=2}^n  c_j$.  This is equivalent to enforcing the column sum $\sum_{i=1}^m x_i =  c_1$.  Choosing the permutation $\boldsymbol \phi$ so that $ r_{\phi_1}-x_{\phi_1}\geq\dotsb\geq  r_{\phi_m}-x_{\phi_m}$, the remaining Gale--Ryser conditions are that
\begin{equation} \sum_{\ell=1}^i ( r_{\phi_\ell}-x_{\phi_\ell}) \leq \sum_{\ell=1}^i \sum_{j\geq 2} {\mathds{1}}\{ c_j\geq \ell\} \quad \text{for all $i=1,\dotsc,m-1$}, \label{e:Bphi} \end{equation}
which implies that
\begin{equation}  \sum_{\ell=1}^i ( r_{\pi_\ell}-x_{\pi_\ell}) \leq \sum_{\ell=1}^i \sum_{j\geq 2} {\mathds{1}}\{ c_j\geq \ell\} \quad \text{for all $i=1,\dotsc,m-1$}, \label{e:Bpi} \end{equation}
since the permutation $\boldsymbol \phi$ makes the left side as large as possible.  Solving \eqref{e:Bpi} for $\sum_{\ell=1}^i x_{\pi_\ell}$ gives the bounds encoded in the $\mathcal{B}_i$ and shows that $\tilde\Omega\supseteq\Omega$.

We cannot use the permutation $\boldsymbol \phi$ in the construction of $Q$ because $\boldsymbol \phi$ depends on $\boldsymbol x$, however, \citet{Chen:Sequential:2005} further prove that \eqref{e:Bphi} and \eqref{e:Bpi} are in fact equivalent, which means we are in the ideal situation where $\Omega=\tilde\Omega$.  (Although they made use of the factorization in Theorem \ref{t:binary}, their proposal distributions were not of the form in \eqref{e:Q}, except in the special case where $\Omega=\{\boldsymbol x\in\{0,1\}^{m\times 1}:\sum_i x_i=c_1\}$.)  Furthermore, \citet{chen2007cit} provides an extension of Theorem \ref{t:binary} for the case where, in addition to the margin constraints, $\Omega^*$ also enforces a fixed pattern of structural zeros for which there is at most one structural zero in each row and column.  This includes the important special case of adjacency matrices of directed graphs; see supplementary material.

\subsection{Combinatorial approximations} \label{s:U}

Here we discuss approximation of $U$ by $\tilde U$.  Define
\[ N_{m,n}(\boldsymbol r,\boldsymbol c) = \bigl|\{\boldsymbol z\in\{0,1\}^{m\times n}:\boldsymbol{R}(\boldsymbol z)=\boldsymbol r, \ \boldsymbol C(\boldsymbol z)=\boldsymbol c\}\bigr| \]
to be the number of $m\times n$ binary matrices with row sums $\boldsymbol r$ and column sums $\boldsymbol c$.  We have
\[ U(\boldsymbol x) = \Prob(\boldsymbol Y^1=\boldsymbol x) = \frac{N_{m,n-1}(\boldsymbol r-\boldsymbol x,\boldsymbol c^{2:n})}{N_{m,n}(\boldsymbol r,\boldsymbol c)} , \]
(where as before, $\boldsymbol Y$ is uniform over $\Omega^*$)
and we desire an approximation of the form
\begin{equation} U(\boldsymbol x) = \frac{N_{m,n-1}(\boldsymbol r-\boldsymbol x,\boldsymbol c^{2:n})}{N_{m,n}(\boldsymbol r,\boldsymbol c)} \approx \gamma \prod_{i=1}^m u_i^{x_i} , \label{e:Nprod} \end{equation}
where $\gamma$ is an irrelevant positive constant.

Temporarily pretending that \eqref{e:Nprod} is accurate for any $\boldsymbol x\in\{0,1\}^{m\times 1}$, we have
\[ u_i \approx \frac{U(\boldsymbol I^i)}{U(\boldsymbol 0)} =  \frac{N_{m,n-1}(\boldsymbol r-\boldsymbol I^i,\boldsymbol c^{2:n})}{N_{m,n-1}(\boldsymbol r,\boldsymbol c^{2:n})} , \]
where $\boldsymbol I^i$ is the $i$th column of the $m\times m$ identity matrix $\boldsymbol I$.  We cannot use this directly, since it is trying to evaluate $U$ outside of $\Omega$, and, furthermore, computationally efficient procedures for evaluating $N_{m,n-1}$ are not available.  Nevertheless, it suggests using
\begin{equation} u_i = \tilde N_{m,n-1}(\boldsymbol r-\boldsymbol I^i,\boldsymbol c^{2:n})/\tilde N_{m,n-1}(\boldsymbol r,\boldsymbol c^{2:n}) \quad \quad (i=1,\dotsc,m)  \label{e:u} \end{equation}
for an approximation $\tilde N$ of $N$ that extends smoothly to invalid margins.

Several asymptotic approximations for $N$ have appeared in the literature and could be used for $\tilde N$.  For example, \citet[Theorem 1]{canfield2008aed} suggest the following, which we write asymmetrically with respect to $\boldsymbol r$ and $\boldsymbol c$ in order to simplify \eqref{e:canfield-u} below:
\[ \begin{gathered}  \tilde N_{m,n}(\boldsymbol r,\boldsymbol c) = \binom{mn}{\sum_{k=1}^n c_k}^{-1}\prod_{i=1}^m \binom{n}{r_i}\prod_{j=1}^n \binom{m}{c_j} \exp\Bigl(-\frac{1}{2}\bigl(1-\mu_{m,n}(\boldsymbol r,\boldsymbol c)\bigr)\bigl(1-\nu_{m,n}(\boldsymbol c)\bigr)\Bigr) , \\
\mu_{m,n}(\boldsymbol r,\boldsymbol c) = \eta_{m,n}(\boldsymbol c)\sum_{i=1}^m \biggl(r_i-\frac{1}{m}\sum_{k=1}^n c_k\biggr)^2 , \quad \quad \nu_{m,n}(\boldsymbol c) = \eta_{m,n}(\boldsymbol c)\sum_{j=1}^n \biggl(c_j-\frac{1}{n}\sum_{k=1}^n c_k\biggr)^2  , \\
\eta_{m,n}(\boldsymbol c) = \frac{mn}{(\sum_{k=1}^n c_k) (mn-\sum_{k=1}^n c_k)} .
\end{gathered} \]
Substituting this into \eqref{e:u} and simplifying gives
\begin{equation} \label{e:canfield-u} u_i = \frac{r_i}{n-r_i}\exp\biggl[\eta_{m,n-1}(\boldsymbol c^{2:n})\bigl(1-\nu_{m,n-1}(\boldsymbol c^{2:n})\bigr)\biggl(\frac{1}{2}-r_i+\frac{1}{m}\sum_{k=2}^n c_k \biggr)\biggr] . \end{equation}
If $r_i=0$ or $r_i=n$, then the value of $X_i$ is determined by $\tilde\Omega$, and any choice of $u_i>0$ gives the same $Q$; we use $u_i=1$ in these cases.
We find that \eqref{e:canfield-u} works well over a large range of margins when $P^*$ is uniform.  It is excellent if the margins are approximately semi-regular, that is, if the row and column sums do not deviate substantially from their respective mean values.

For certain pathological cases with wildly varying margins, such as those in \citet{bezakova2006nes}, \eqref{e:canfield-u} does not work well.  However, if the margins are such that the resulting matrices have a very low density of ones, even if the margins are highly irregular, then good performance can be obtained by instead using the asymptotic approximation of $N$ from \citet[Theorem 1.3]{greenhill2006aes}.  Details are provided in the supplementary material.  In fact, for the specific pathological cases in \citet{bezakova2006nes} using this alternative approximation gives $Q^*\equiv P^*$ in the uniform case.  None of the computationally efficient combinatorial approximations that we have found in literature work well when the margins are both highly irregular and lead to a moderate density of ones, but we are hopeful that advances in asymptotic enumeration techniques will eventually lead to approximations that work well in almost all cases.

\citet{Chen:Sequential:2005} observed that combinatorial approximations could be used to find a good choice of $\boldsymbol u$ and they mentioned an early asymptotic approximation from \citet{Oneil:Asymptotics:1969}, which was explored further by \citet{blanchet2006isa} in an asymptotic analysis of the algorithm.  The examples in \citet{Chen:Sequential:2005}, however, use $u_i=r_i/(n-r_i)$, which is motivated by considering only the row margin constraints.  Although there are several substantial differences between their proposal distribution and ours for the special case of the uniform distribution of $\Omega^*$, we suspect that much of the improved performance of our algorithm results from using more accurate combinatorial approximations.

In the next section we use $V$ and $\tilde V$ to account for the effects of $\boldsymbol w$, including the effects of zeros in $\boldsymbol w$.  Since these zeros affect the size of the support of $P^*$, an alternative, perhaps more natural approach is to allow $U$ and $\tilde U$ to capture the effects of zeros in $\boldsymbol w$.  The supplementary material contains more details.

\subsection{Non-uniform weighting} \label{s:V}

Here we discuss approximation of $V$ by $\tilde V$.  To develop an approximation, we will ignore the column margins and consider only the row margins.  This is similar to the approach used by \citet{Chen:Sequential:2005} to develop combinatorial approximations for the uniform case.  Let $\boldsymbol B\in\{0,1\}^{m\times n}$ be a matrix of independent Bernoulli$(1/2)$ random variables, so that
\begin{align} V(\boldsymbol x) & = \Exp\Bigl(\prod_{ij} w_{ij}^{Y_{ij}} \Bigl| \boldsymbol Y^1=\boldsymbol x\Bigr) = \Exp\Bigl(\prod_{ij} w_{ij}^{B_{ij}} \Bigl| \boldsymbol B^1=\boldsymbol x, \boldsymbol R(\boldsymbol B)=\boldsymbol r, \boldsymbol C(\boldsymbol B)=\boldsymbol c\Bigr) \notag
\\
&
\approx \Exp\Bigl(\prod_{ij} w_{ij}^{B_{ij}} \Bigl| \boldsymbol B^1=\boldsymbol x, \boldsymbol R(\boldsymbol B)=\boldsymbol r\Bigr) = \prod_{i=1}^m \Exp\Bigl(\prod_{j=1}^n w_{ij}^{B_{ij}} \Bigl| B_{i1}=x_i, R_i(\boldsymbol B)=r_i\Bigr) \notag
\\
&
\propto \prod_{i=1}^m v_i^{x_i} ,
\label{e:tildeV} \end{align}
where
\begin{align} v_i & = \frac{\Exp\bigl(\prod_{j=1}^n w_{ij}^{B_{ij}} \bigl| B_{i1}=1, R_i(\boldsymbol B)=r_i\bigr)}{\Exp\bigl(\prod_{j=1}^n w_{ij}^{B_{ij}} \bigl| B_{i1}=0, R_i(\boldsymbol B)=r_i\bigr)}  \quad \quad (r_i=1,\dotsc,n-1; \ i=1,\dotsc,m) \notag
\\
& = \frac{w_{i1}\binom{n-1}{r_i-1}^{-1}\sum_{\boldsymbol b\in\{0,1\}^{n-1}} {\mathds{1}}\{{\textstyle \sum_{j=1}^{n-1} b_j=r_i-1}\}\prod_{j=2}^n w_{ij}^{b_{j-1}}}{\binom{n-1}{r_i}^{-1}\sum_{\boldsymbol b\in\{0,1\}^{n-1}} {\mathds{1}}\{{\textstyle \sum_{j=1}^{n-1} b_j=r_i}\}\prod_{j=2}^n  w_{ij}^{b_{j-1}}} .
\label{e:v1}
\end{align}
In the supplementary material we describe how to compute all possible $\boldsymbol v$ for all columns using $O(nd)$ operations (where $d=\sum_i r_i=\sum_j c_j$) in a one-time preprocessing step that can be done prior to sampling.  As with $u_i$, we always define $v_i=1$ if $r_i=0$ or $r_i=n$.  For cases where $\boldsymbol w$ has zeros, we can sometimes have a zero in the denominator of \eqref{e:v1} for $r_i>0$.  This happens when fewer than $r_i$ of the $n-1$ remaining weights in the row are nonzero.  Consequently, we need to force $X_i=1$, which we do by setting the corresponding $\mathcal A_{\ell}=\{1\}$ in Section \ref{s:margin}.

An important observation that we have thus far neglected is that many different choices of $\boldsymbol w$ give rise to the same $P^*$.  Define
\[ \Lambda(\boldsymbol w) = \{\boldsymbol t \in [0,\infty)^{m\times n} : \boldsymbol t = \boldsymbol \alpha\boldsymbol \beta^{\text{\sc t}}\circ \boldsymbol w,\ \boldsymbol \alpha\in(0,\infty)^{m\times 1}, \ \boldsymbol \beta\in(0,\infty)^{n\times 1}\} , \]
where ${\text{\sc t}}$ denotes transpose and $\boldsymbol y\circ \boldsymbol z$ is the Hadamard product, that is, element-wise multiplication of matrices of the same size, defined by $(\boldsymbol y\circ \boldsymbol z)_{ij}=y_{ij}z_{ij}$.
Then for every $\boldsymbol t\in \Lambda(\boldsymbol w)$ it is straightforward to verify that the $P^*$ defined with the weight matrix $\boldsymbol w$ and the $P^*$ defined with the weight matrix $\boldsymbol t$ are identical.  Similarly, the two versions of $V$ differ only by an inconsequential constant of proportionality.  Unfortunately, our approximation $\tilde V(\boldsymbol x)\propto\prod_i v_i^{x_i}$ defined above does not share this invariance.  Consequently, proposal distributions constructed with $\boldsymbol w$ and $\boldsymbol t$, respectively, could differ, even though the target distribution does not differ.
We find this unappealing and remedy it in a preprocessing step prior to the construction of $Q^*$ by first transforming $\boldsymbol w$ into an equivalent, canonical $\bar{\boldsymbol w}\in\Lambda(\boldsymbol w)$.  In particular, $\bar{\boldsymbol w}$ is the unique element of $\Lambda(\boldsymbol w)$ with the property that its average nonzero entry over any row or column is one, namely,
\begin{equation} \sum_{i=1}^m \bar w_{ij} = \sum_{i=1}^m {\mathds{1}}\{\bar w_{ij} > 0\} ,\quad  \sum_{j=1}^n \bar w_{ij} = \sum_{j=1}^n {\mathds{1}}\{\bar w_{ij} > 0\}  \quad (i=1,\dotsc,m;\ j=1,\dotsc,n) . \label{e:barwH} \end{equation}
The solution to \eqref{e:barwH} over $\Lambda(\boldsymbol w)$ exists, is unique, and is easy to find numerically \citep{rothblum1989scalings}; see supplementary material for details and more discussion.  In the examples below, we always define $\boldsymbol v$ in terms of $\bar{\boldsymbol w}$, not $\boldsymbol w$.  Not only does this ensure that $Q^*$ has the same invariance property as $P^*$, but we find that performance of the algorithm tends to improve.

If we know that $P^*$ is uniform over $\Omega^*$, for example, if $\boldsymbol w\equiv 1$ or $\bar{\boldsymbol w}\equiv 1$, then $V$ is constant over $\Omega$, and we can ignore $\boldsymbol v$ in the construction of $Q$.

\section{Importance sampling}

\subsection{Algorithm summary} \label{s:overview}

\begin{enumerate}
\item  \label{l:pre} Preprocessing: Compute $\bar{\boldsymbol w}$ from \eqref{e:barwH} and precompute all possible $\boldsymbol v$ using \eqref{e:v1} with $\bar{\boldsymbol w}$ in place of $\boldsymbol w$; see supplementary material for details.  Compute $\boldsymbol\pi$, $\boldsymbol{\mathcal{A}}$, and $\boldsymbol{\mathcal{B}}$ according to Theorem \ref{t:binary} and $\boldsymbol{u}$ according to \eqref{e:canfield-u} for the first column.
\item \label{l:loop} Generating a single observation, $\boldsymbol{\tilde Z}$, from $Q^*$:  The matrix $\boldsymbol{\tilde Z}$ is generated column-by-column as follows. Set $q=1$. Sequentially, for each column:
\begin{enumerate}
\item \label{2a} For the current $m,n,\boldsymbol r,\boldsymbol c$, compute $\boldsymbol\pi,\boldsymbol{\mathcal{A}},\boldsymbol{\mathcal{B}},\boldsymbol u$ as above.  After the preprocessing, updating these quantities based on the previously sampled column requires $O(m)$ operations.
\item Use Theorem \ref{t:DP} to compute the Markov chain representation for $Q$.  For the $j$th column, this takes $O(mc_j)$ operations.  If $Q\equiv 0$, then $\boldsymbol{\tilde Z}$ will not be in the support of $P^*$; assign a final importance weight of zero and go to step \ref{l:done}.
\item \label{l:makeX} Generate a random observation $\boldsymbol X$ from $Q$ and evaluate $Q(\boldsymbol X)$.  This takes $O(m)$ operations.
\item Assign the current column of $\boldsymbol{\tilde Z}$ to be $\boldsymbol X$.  Update $q \leftarrow qQ(\boldsymbol X)$, $\boldsymbol c \leftarrow \boldsymbol c^{2:n}$, $n\leftarrow n-1$, $\boldsymbol r\leftarrow \boldsymbol r-\boldsymbol X$.  If $n >0$, continue looping over columns; go to step \ref{2a}.  If $n=0$, the final matrix is $\boldsymbol{\tilde Z}$, and we have $Q^*(\boldsymbol{\tilde Z})=q$; go to step \ref{l:done}.
\end{enumerate}
\item \label{l:done} To generate additional independent observations from $Q^*$, reset all variables to their original values after step \ref{l:pre} and repeat step \ref{l:loop}.
\end{enumerate}
The same algorithm can be used to evaluate $Q^*(\boldsymbol z)$ for any $\boldsymbol z\in\Omega^*$.  Simply assign $\boldsymbol X$ to be the current column of $\boldsymbol z$ in step \ref{l:makeX}, instead of sampling a new column.  (The algorithm can be applied for any ordering of the columns, and $Q^*$ will depend on the chosen ordering.  The supplementary material describes the heuristics that we use to choose a column ordering.)

\subsection{Monte Carlo approximation and diagnostics}

Let $\boldsymbol Z$ have distribution $P^*$, let $\boldsymbol Z_1,\dotsc,\boldsymbol Z_T$ be random sample from $Q^*$ generated as above, and let $h$ be a function over $\Omega^*$.  Define the unnormalized importance weights
\[ f(\boldsymbol z) = \frac{\kappa P^*(\boldsymbol z)}{Q^*(\boldsymbol z)} = \frac{\prod_{ij} w_{ij}^{z_{ij}}{\mathds{1}}\{\boldsymbol z\in\Omega^*\}}{Q^*(\boldsymbol z)}  \quad \quad (\boldsymbol z\in\{0,1\}^{m\times n}) , \]
which we can efficiently evaluate for any $\boldsymbol z$ as described above.  In the formula for $f(\boldsymbol z)$ it is important that we use $\boldsymbol w$, not $\bar{\boldsymbol  w}$, particularly if we are approximating $\kappa$.
We can approximate $\kappa$ and $\mu=\Exp(h(\boldsymbol Z))$ via importance sampling in the usual way, namely,
\begin{equation} \hat\mu_T = \frac{\sum_{t=1}^T f(\boldsymbol Z_t) h(\boldsymbol Z_t)}{\sum_{t=1}^T f(\boldsymbol Z_t)} \to \mu , \quad\quad \hat\kappa_T = \frac{1}{T}\sum_{t=1}^T f(\boldsymbol Z_t) \to \kappa \quad \quad (T\to\infty) .  \label{e:kappahat} \end{equation}
Besides being consistent, $\hat\kappa_T$ and $\hat\kappa_T\hat\mu_T$ are also unbiased approximations of $\kappa$ and $\kappa\mu$, respectively.  See \citet{Liu:Monte:2001} for details about importance sampling.  See \citet{harrison2012conservative} for modifications when \eqref{e:kappahat} is used to approximate a p-value.

In this context, importance sampling algorithms are usually evaluated empirically by diagnostics related to the variability of the importance weights.  The less variable the importance weights, the better the algorithm is judged to be performing.  For the numerical illustrations below, we report
\[ \hat\text{\small\sc cv}_T^2 = \frac{1}{(T-1)\hat\kappa_T^2}\sum_{t=1}^T (f(\boldsymbol Z_t)-\hat\kappa_T)^2 , \quad \quad \hat \Delta_T = \frac{\max_{t=1,\dotsc,T} f(\boldsymbol Z_t)}{\min_{t=1,\dotsc,T} f(\boldsymbol Z_t)}-1. \]
As $T\to\infty$, the approximate squared coefficient of variation, $\hat\text{\small\sc cv}_T^2$, converges to the true squared coefficient of variation, $\text{\small\sc cv}^2=\Var(f(\boldsymbol Z_1))/\Exp(f(\boldsymbol Z_1))^2=\Var(P^*(\boldsymbol Z_1)/Q^*(\boldsymbol Z_1))$, where $\Var$ denotes variance.  $\tilde T = T/(1+\text{\small\sc cv}^2)$ has been suggested as a rough diagnostic for effective sample size, meaning that a sample size of $T$ from $Q^*$ behaves roughly like a sample size of $\tilde T$ from $P^*$ for the purposes of Monte Carlo approximating $\mu$ for well-behaved functions $h$ \citep{Kong:Sequential:1994,Liu:Monte:2001}.  For many but not all examples we find $\hat\text{\small\sc cv}^2_T < 1$, suggesting that $Q^*$ is appropriate for efficient importance sampling.  The relative range of importance weights reported by $\hat \Delta_T$ is an especially stringent diagnostic.  For nearly constant margins and $P^*$ close to uniform, we often find $\hat \Delta_T\approx 0$, suggesting that $Q^*$ is an excellent approximation of $P^*$; see Table \ref{t:1}.

\section{Numerical illustrations} \label{s:sim}

We experiment with four different classes of weights based on a canonical matrix $\boldsymbol y$ whose entries are independently sampled from the uniform$(0,1)$ distribution: (I) $w_{ij}= 1$, which is the uniform distribution over $\Omega^*$, (II) $w_{ij}=y_{ij}+1$, (III) $w_{ij}=y_{ij}$, and (IV) $w_{ij}=-{\mathds{1}}(y_{ij}<0.99)\log(y_{ij})$, for all $i,j$. The specific entries of $\boldsymbol y$ for different sized matrices are in the supplementary material.  The resulting $P^*$ is increasingly non-uniform in each of the latter three cases and has $1\%$ structural zeros in case (IV).  Recall that each $\boldsymbol w$ corresponds to a family of weights of the form $\boldsymbol \alpha\boldsymbol \beta^{\text{\sc t}}\circ \boldsymbol w$ that give the same $P^*$ and $Q^*$; see Section \ref{s:V}.  In all cases we report results with $T=1000$.

We begin with $500\times 500$ $r_1$-regular matrices, i.e., $r_i=c_j=r_1$ for all $i,j=1,\dotsc,500$.  Results are summarized in Table \ref{t:1} for $r_1=1,2,4,8,\dotsc,256$.  The diagnostics are striking, especially in the uniform case, for which the importance weights are essentially constant.
Performance degrades slightly as $P^*$ becomes strongly non-uniform, but in all cases the estimated $\text{\small\sc cv}^2$ is less than one.  Low variability in importance weights corresponds to high precision in estimates of $\kappa$.  For example, in the uniform case, where $\kappa=|\Omega^*|$, for $r_1=256$ we obtain $\hat\kappa = (1.478301\pm0.000044){\times}10^{73781}$, where the errors are approximate standard errors estimated from the same importance samples,  and for $r_1=2$ we obtain $\hat\kappa_T = (2.27653\pm0.00017){\times}10^{2266}$, the latter of which is close to the true value of $\kappa=2.27658\dotsc{\times}10^{2266}$; see supplementary material.  To our knowledge the exact value of $\kappa$ in these examples can only be efficiently computed for the special case of the uniform distribution over either $1$-regular or $2$-regular matrices \citep{anand1966combinatorial}.  Sampling from the uniform distribution over $1$-regular matrices is trivial, $\kappa=m!$,  and there is no need to use our algorithm, although it is comforting that $Q^*=P^*$ in this case.

We remark that the distributions corresponding to different weight classes in Table \ref{t:1} are almost singular with respect to each other.  For example, in the $1$-regular case, if we use the $Q^*$ for weight class I as a proposal distribution for the $P^*$ corresponding to one of the other weight classes, then we obtain, for weight class II, $\hat\Delta_T=2{\times}10^{13}$ and $\hat\text{\small\sc cv}_T^2=8{\times}10^2$, for weight class III, $\hat\Delta_T=2{\times}10^{67}$ and $\hat\text{\small\sc cv}_T^2=1{\times}10^3$, and for weight class IV, only $6$ of the $1000$ observations from $Q^*$ were even in the support of $P^*$, owing to the structural zeros.  Results are similar for other combinations and become even more extreme as $r_1$ increases.

\begin{table}[h]
\centering
\caption{$500\times 500$ $r_1$-regular matrices}
{\footnotesize
\begin{tabular}{ccccccccccccc}
& & \multicolumn{2}{c}{uniform} & & \multicolumn{2}{c}{$\boldsymbol w$ class II} & & \multicolumn{2}{c}{$\boldsymbol w$ class III} & & \multicolumn{2}{c}{$\boldsymbol w$ class IV} \\ [1ex]
$r_1$$\vphantom{\hat{\bigl(}}$ & & $\hat\Delta_T$ & $\hat\text{\small\sc cv}^2_T$ & & $\hat\Delta_T$ & $\hat\text{\small\sc cv}^2_T$ & & $\hat\Delta_T$ & $\hat\text{\small\sc cv}^2_T$ & & $\hat\Delta_T$ & $\hat\text{\small\sc cv}^2_T$ \\ [1ex]
$1$  & &  $0$ &  $0$ & &  $2{\times}10^{-1}$ &  $5{\times}10^{-4}$ & &  $4{\times}10^{0}$ &  $4{\times}10^{-2}$ & &  $5{\times}10^{1}$ &  $3{\times}10^{-1}$ \\
 $2$  & &  $4{\times}10^{-2}$ &  $5{\times}10^{-6}$ & &  $2{\times}10^{-1}$ &  $4{\times}10^{-4}$ & &  $6{\times}10^{0}$ &  $4{\times}10^{-2}$ & &  $8{\times}10^{1}$ &  $2{\times}10^{-1}$ \\
 $4$  & &  $1{\times}10^{-2}$ &  $1{\times}10^{-6}$ & &  $1{\times}10^{-1}$ &  $4{\times}10^{-4}$ & &  $5{\times}10^{0}$ &  $3{\times}10^{-2}$ & &  $2{\times}10^{2}$ &  $2{\times}10^{-1}$ \\
 $8$  & &  $2{\times}10^{-2}$ &  $1{\times}10^{-6}$ & &  $2{\times}10^{-1}$ &  $3{\times}10^{-4}$ & &  $3{\times}10^{0}$ &  $3{\times}10^{-2}$ & &  $4{\times}10^{1}$ &  $2{\times}10^{-1}$ \\
 $16$  & &  $1{\times}10^{-2}$ &  $1{\times}10^{-6}$ & &  $2{\times}10^{-1}$ &  $3{\times}10^{-4}$ & &  $3{\times}10^{0}$ &  $3{\times}10^{-2}$ & &  $4{\times}10^{1}$ &  $1{\times}10^{-1}$ \\
 $32$  & &  $8{\times}10^{-3}$ &  $8{\times}10^{-7}$ & &  $1{\times}10^{-1}$ &  $2{\times}10^{-4}$ & &  $2{\times}10^{0}$ &  $2{\times}10^{-2}$ & &  $1{\times}10^{1}$ &  $1{\times}10^{-1}$ \\
 $64$  & &  $9{\times}10^{-3}$ &  $9{\times}10^{-7}$ & &  $1{\times}10^{-1}$ &  $2{\times}10^{-4}$ & &  $3{\times}10^{0}$ &  $2{\times}10^{-2}$ & &  $2{\times}10^{1}$ &  $9{\times}10^{-2}$ \\
 $128$  & &  $1{\times}10^{-2}$ &  $9{\times}10^{-7}$ & &  $1{\times}10^{-1}$ &  $9{\times}10^{-5}$ & &  $1{\times}10^{0}$ &  $1{\times}10^{-2}$ & &  $5{\times}10^{0}$ &  $5{\times}10^{-2}$ \\
 $256$  & &  $9{\times}10^{-3}$ &  $9{\times}10^{-7}$ & &  $5{\times}10^{-2}$ &  $5{\times}10^{-5}$ & &  $1{\times}10^{0}$ &  $1{\times}10^{-2}$ & &  $9{\times}10^{0}$ &  $7{\times}10^{-2}$
\end{tabular}
}
\label{t:1}
\end{table}

For the special case of $1$-regular matrices, corresponding to the first row in Table \ref{t:1}, $\kappa$ is the permanent of $\boldsymbol w$ and various generalizations of the permanent correspond to expectations under $P^*$.  The current state-of-the-art algorithm for approximating permanents and $\alpha$-permanents of general matrices, see \eqref{e:perm} above, seems to be the importance sampling algorithm of \citet{kou2009approximating}, which has about the same computational complexity as our algorithm.  For the case $\alpha=1$, their algorithm is nearly identical to ours, the main differences being the choice of column order and our use of $\boldsymbol{\bar w}$, which seems to give our algorithm slightly better performance.  Their algorithm is generally preferable for $\alpha\neq 1$, since it is tailored to the specific choice of $\alpha$, although in many cases performance is comparable.  The supplementary materials have numerical comparisons for each of the $\boldsymbol w$ used in Table \ref{t:1} and for all of the examples in \citet{kou2009approximating}, which include cases with $\alpha=1$ and $\alpha=1/2$.  It is interesting that in many cases our generic approach is competitive with specialized software.

In Table \ref{t:2} we repeat the simulations of Table \ref{t:1} for $50\times 100$ irregular matrices with margins $\boldsymbol r=k\tilde{\boldsymbol  r}$ and $\boldsymbol c=k\tilde{\boldsymbol  c}$, for the cases $k=1,\dotsc,4$, where ${\tilde{\boldsymbol  r}}^{\text{\sc t}}= (24^1$, $22^2$, $17^4$, $13^3$, $12^2$, $11^3$, $10^2$, $9^3$, $8^6$, $7^1$, $6^4$, $5^4$, $4^5$, $3^6$, $2^4)$ and ${\tilde{\boldsymbol  c}}= (12^2$, $10^2$, $9^5$, $8^4$, $7^6$, $6^{11}$, $5^{10}$, $4^{18}$, $3^9$, $2^{13}$, $1^{20})$ using $i^j$ to denote $j$ copies of $i$.  Performance degrades in the irregular case as the matrices become more dense.  In most, but not all cases, the diagnostics suggest $Q^*$ could be used for efficient importance sampling.

\begin{table}[h]
\centering
\caption{$50\times 100$ irregular matrices with $\boldsymbol r=k\tilde{\boldsymbol  r}$, $\boldsymbol c=k\tilde{\boldsymbol  c}$}
{\footnotesize
\begin{tabular}{ccccccccccccc}
& & \multicolumn{2}{c}{uniform} & & \multicolumn{2}{c}{$\boldsymbol w$ class II} & & \multicolumn{2}{c}{$\boldsymbol w$ class III} & & \multicolumn{2}{c}{$\boldsymbol w$ class IV} \\ [1ex]
$k$$\vphantom{\hat{\bigl(}}$ & & $\hat\Delta_T$ & $\hat\text{\small\sc cv}^2_T$ & & $\hat\Delta_T$ & $\hat\text{\small\sc cv}^2_T$ & & $\hat\Delta_T$ & $\hat\text{\small\sc cv}^2_T$ & & $\hat\Delta_T$ & $\hat\text{\small\sc cv}^2_T$ \\ [1ex]
 $1$
  & &  $4{\times}10^{-1}$ &  $1{\times}10^{-3}$
  & &  $3{\times}10^{0\phantom{}}$ &  $5{\times}10^{-2}$
  & &  $8{\times}10^{1\phantom{}}$ &  $5{\times}10^{-1}$
  & &  $5{\times}10^{3\phantom{1}}$ &  $3{\times}10^{0\phantom{}}$
 \\
 $2$
  & &  $3{\times}10^{0\phantom{-}}$ &  $3{\times}10^{-2}$
  & &  $7{\times}10^{0\phantom{}}$ &  $1{\times}10^{-1}$
  & &  $7{\times}10^{2\phantom{}}$ &  $2{\times}10^{0\phantom{-}}$
  & &  $6{\times}10^{4\phantom{1}}$ &  $7{\times}10^{0\phantom{}}$
 \\
 $3$
  & &  $2{\times}10^{2\phantom{-}}$ &  $7{\times}10^{-1}$
  & &  $2{\times}10^{2\phantom{}}$ &  $6{\times}10^{-1}$
  & &  $2{\times}10^{4\phantom{}}$ &  $6{\times}10^{0\phantom{-}}$
  & &  $3{\times}10^{6\phantom{1}}$ &  $4{\times}10^{1\phantom{}}$
 \\
 $4$
  & &  $3{\times}10^{6\phantom{-}}$ &  $3{\times}10^{1\phantom{-}}$
  & &  $3{\times}10^{6\phantom{}}$ &  $2{\times}10^{1\phantom{-}}$
  & &  $4{\times}10^{9\phantom{}}$ &  $2{\times}10^{2\phantom{-}}$
  & &  $2{\times}10^{13}$ &  $8{\times}10^{2\phantom{}}$
\end{tabular}
}
\label{t:2}
\end{table}

For the special case of the uniform distribution, corresponding to the far left category of weights in Tables \ref{t:1} and \ref{t:2}, the sequential importance sampling algorithm of \citet{Chen:Sequential:2005}, as implemented in the publicly available R package {\tt networksis} \citep{admiraal2008networksis}, appears to be the current state-of-the-art algorithm for practical Monte Carlo approximation.  Our algorithm is a substantial improvement, especially for dense or irregular margins.  Using {\tt networksis} gives $\hat\Delta_T=(1{\times}10^1,1{\times}10^2,2{\times}10^5,1{\times}10^{11})$ and $\hat\text{\small\sc cv}^2_T=(2{\times}10^{-1},6{\times}10^{-1},1{\times}10^1,4{\times}10^2)$ for the first pair of columns in Table \ref{t:2}.  The {\tt networksis} implementation is several orders of magnitude slower than our implementation, and is too slow for most of the examples in Table \ref{t:1}.

redpajama_set_name: RedPajamaArXiv


text:
Az East Ham a londoni metró egyik állomása a 3-as és 4-es zóna határán, a District line és a Hammersmith & City line érinti.

Története
Az állomást 1858-ban adták át a London, Tilbury and Southend Railway részeként, melyet 1902-től a District Railway vonatai szolgálták ki. A Metropolitan Railway 1936-tól használta az állomást, majd 1990-ben átnevezték a vonalat és a Hammersmith & City line nevet kapta.

Forgalom

Átszállási kapcsolatok

Fordítás

Források

Londoni metróállomások

redpajama_set_name: RedPajamaWikipedia


You may have spotted a few green beauty posts lately on my blog. Even though I have no preference for natural over synthetics, I do like the idea of using natural, organic ingredients. That's why I was so thrilled when I encountered Susanne Kaufmann: an Austrian brand with an intriguing history. The brand actually started as a spa in Hotel Post Bezau. Together with a team of pharmacists and dermatologists, Susanne started to develop products for the applications in this spa from precious natural ingredients. I honestly wish I could visit their spa, it looks totally amazing! If I only had the time for a little get-away… Now, let's move on to the review of two products: the Shower Line M & Refreshing Tea.
The Shower/Shampoo Line M (250 ml) was specifically designed for men. According to the brand, it's a shampoo, conditioner, shower gel and skin nourishment all in one. The first thing that I have to mention is the packaging – the amber coloured glass bottle and black logo make this product look incredibly chique in my shower. The addition of a pump makes it easy to dose the product rightfully – no fussing with the pump and no spillage!
When I tried this the first time, I was struck by the gorgeous scent: a deep, yet clear musk-y scent that reminds me of a fresh, natural cologne. In this case, I'm confident to say that most men would enjoy this scent. As this is "specifically designed" for men, there's also a unisex variant available. The product itself is quite low-foaming and has a light texture to it. I know some of you dislike low-foaming products, but I actually enjoy using them more as they leave my skin soft & supple instead of incredibly dry. After using this, I feel no need to slather on a rich body lotion to soothe my dry patches.
I also tried it as a shampoo and conditioner for my hair and it worked fine! It left my hair feeling clean, moisturized and not dry at all. However, I am literally restriced to using specific hair care products for my grey hair as I can't take care of it with just one product. For those who don't have any concerns, this works brilliantly.
As the bottle itself is made of glass, you really have to be careful with it! Even though it's aesthetically pleasing, it might be difficult for some people (read: me) to not drop and break it. On the other hand, it's much more environmentally friendly than a plastic bottle!
Ingredients: acqua (water), sodium C14-C16 olefin sulfonate, cocamidopropyl betaine, polyglyceryl-3 caprate, glycerin, pentylene glycol, rosa canina (rose hip) extract/oil, hamamelis virginia (witch hazel) distillate/extract, hydroxypropyl guar hydroxypropyltrimonium chloride, parfum (fragrance), sericin, sodium benzoate, tocopherol, alcohol plant, linalool, limonene, geraniol, eugenol, coumarin, citronellol.
How adorable is the little birdy on the packaging?! In fact, the little birdy is their logo – too cute, really. I found out that Susanne's husband's last name is Vogel (German for bird) and she wanted to have a connection between her family and the brand Susanne Kaufmann. Isn't that beautiful?
The Refreshing Tea contains 15 mountain tea bags and it's recommended to use 1 tea bag for 1 litre can and let it steep for about 6 minutes. According to Susanne, the refreshing tea invigorates the senses, freshens the breath and brings new strength and energy. I brew the tea with the instructions in mind and took my first sip – I noticed the lemon and marigold strongly. I was expecting some kind of mint, as it's included in the ingredients, but it isn't very noticeable. I really like it, even though it's not as strong as I'm used to with other teas. It tastes wonderfully and actually does freshen your breath! If you dislike any kind of mint and still want a refreshing tea, then this is the one for you.
Ingredients: marigold leaves, honeybush, lemongrass, verbena, crisped mint, peppermint.
I'm positively surprised by the quality of the products! Susanne Kaufmann really does her best to only include the highest quality of ingredients and make everything smell wonderful. I discovered that skin care is their strongest focus, so I might have to give that a try as well. The Shower/Shampoo Line M (250 ml) and Refreshing Tea (19,5 g) are both available in the Susanne Kaufmann webshop and cost 48,00 euro and 16,00 euro respectively.
Have you tried any organic beauty products recently?
PS: Don't forget to enter my current giveaway!


text:
HomeNablus
Jewish Settler Kills Palestinian Teenager near Qalqilya
A Palestinian teenager was killed Saturday night by a Jewish settler near the settlement of Kedumim, illegally built on Palestinian land in the occupied West Bank province of Qalqilia, the official news agency WAFA reported. […]
Ten Palestinians Injured as Israeli Occupation Forces Storm Nablus
Israeli forces stormed the West Bank city of Nablus on Thursday morning, injuring ten Palestinians and detaining five others, The New Arab reported. The Israeli forces surrounded the city before storming it and targeting resistance […]
PA Suppresses Protest Calling for Release of Political Prisoners in Nablus (VIDEO)
The Palestinian Authority security forces on Tuesday suppressed a protest in Nablus to demand the release of political prisoners held by the PA, Quds Press news site reported. Dozens of masked members of the PA […]
BREAKING: Palestinian Youth Injured in Nablus Dies of His Wounds
A Palestinian youth, identified as Ahmad Abu Junaid, died on Wednesday of wounds sustained during an Israeli army raid of the Balata refugee camp in Nablus, the Ministry of Health announced. The Ministry said in […]
Israeli Forces Shoot, Critically Injure Young Palestinian Man in Nablus (VIDEOS)
Israeli occupation forces shot and critically injured a young Palestinian man on Wednesday, in the Balata refugee camp, in the occupied West Bank city of Nablus, the official news agency WAFA reported. Ahmad Jibril, the […]
Israeli Forces Kill Palestinian Child near Nablus
A Palestinian child has been killed by Israeli forces during a military raid in the occupied West Bank city of Nablus, according to the Palestinian Ministry of Health. Palestinian official news agency WAFA reported that […]
Israeli Forces Shoot, Injure Young Palestinian Man near Nablus
A young Palestinian man was shot and injured on Friday, during an Israeli army raid in the village of Osrin, south of Nablus, the official Palestinian news agency WAFA reported. Ahmad Jibril, the head of […]
Palestinian Toddler Injured in Settler Attack near Nablus
An 18-month-old female Palestinian toddler was injured Friday evening in an attack by Jewish settlers on her family's car in the village of Usarin, near the occupied West Bank city Nablus, the official Palestinian news […]
Palestinians Call on FIFA to Condemn Israeli Killing of Palestinian Footballer
Palestinian Authority Prime Minister Mohammad Shtayyeh has called on the Federation Internationale de Football Association (FIFA) to condemn the killing of Ahmad Atef Mustafa Daraghma, a Palestinian soccer player, who was killed by Israeli forces […]
Israeli Occupation Forces Kill Palestinian Football Player in Nablus
Israeli forces killed a football player and wounded at least 24 Palestinians during a military raid in the northern West Bank city of Nablus. Ahmad Atef Mustafa Daraghma, 23, a resident of the northern West […]
Lions' Den is Not a Fleeting Phenomenon: On Palestine's Looming Armed Revolt
Just when Israel, and even some Palestinians, began talking about the Lions' Den phenomenon in the past tense, a large number of fighters belonging to the newly-formed Palestinian group marched in the city of Nablus. […]
Jewish Settlers Set Fire to Palestinian Nursery, Shop in West Bank
Jewish settlers set fire to a Palestinian nursery and shop on Monday morning in the village of Al-Sawiya, south of Nablus in the northern West Bank, the official Palestinian news agency WAFA reported. According to […]
Lions' Den Say They Successfully Fought off Israeli Raid in Nablus
The Lions' Den Palestinian armed group said in a statement on Wednesday that its fighters successfully fought off an Israeli raid on the old city of Nablus, earlier in the morning, The New Arab reported. […]
Jewish Settlers Attack Palestinian School near Nablus
Jewish settlers attacked a Palestinian high school near the occupied West Bank city of Nablus on Wednesday, the official Palestinian news agency WAFA reported. Ghassan Daghlas, who monitors settlements activities in the northern West Bank, […]
WATCH: Three Palestinians Injured during Israeli Military Raid in Nablus
A Palestinian was shot with a live bullet on Wednesday, during an Israeli military raid in the occupied West Bank city of Nablus, the official Palestinian news agency WAFA reported. Ahmad Jibril, the head of […]
UN Special Coordinator 'Horrified' by Israeli Execution of Palestinian
Tor Wennesland, the United Nations Special Coordinator for the Middle East Peace Process, said on Friday that he was horrified by the scene of an Israeli soldier killing a Palestinian in the town of Hawwara, […]
WATCH: Israeli Forces Kill Young Palestinian Man near Nablus
Israeli forces shot and killed a young Palestinian man near Nablus, the official Palestinian news agency WAFA reported. The Palestinian Ministry of Health announced in a brief statement that a Palestinian man, later identified as […]
WATCH: Israeli Forces Shoot, Injure Young Palestinian Man in Nablus
Israeli forces Friday morning shot a Palestinian man and detained another in the occupied West Bank city of Nablus, the official Palestinian news agency WAFA reported. WAFA correspondent said that a sizable Israeli army force […]
Palestinian Man Succumbs to Wounds Sustained by Israeli Gunfire in July (VIDEO)
A young Palestinian man succumbed on Wednesday night to the wounds he sustained from Israeli gunfire last July, the official Palestinian news agency WAFA reported. The Palestinian Health Ministry said that Mohammad Herzallah, 30, died […]
Young Palestinian Man Succumbs to Wounds Sustained by Israeli Gunfire in Nablus
A young Palestinian man succumbed Wednesday evening to critical wounds he sustained after being shot by Israeli forces in Nablus, the official Palestinian news agency WAFA reported. The Palestinian Ministry of Health said in a […]
Israeli Occupation Forces Kill Palestinian Child in Nablus (VIDEO)
A Palestinian child was killed Wednesday overnight by Israeli gunfire, during a military raid in the occupied West Bank city of Nablus, the Palestinian Ministry of Health reported. The Ministry said in a brief statement […]
Israeli Occupation Forces Arrest Two Palestinian Schoolchildren
Israeli occupation forces on Tuesday arrested two Palestinian schoolchildren near the occupied West Bank City of Nablus, the official Palestinian news agency WAFA reported. According to WAFA, two Palestinian minors, identified as Yousef Odah and […]
Israeli Forces Kill Palestinian Teenager near Nablus
A Palestinian teenager was killed on Wednesday during an Israeli military raid near Nablus, the official Palestinian news agency WAFA reported. According to the Palestinian Ministry of Health, 15-year-old Mahdi Hashash was shot during a […]
Israeli Occupation Forces Launch Massive Arrest Campaign in Nablus
Israeli occupation forces launched a campaign of raids and searches in the occupied West Bank on Tuesday, arresting a number of Palestinians. Israeli forces also stormed the city of Nablus and broke into several homes. […]
Israeli Forces Raid Palestinian School, Intimidate Students, Teachers
Israeli occupation forces on Sunday raided a school serving the villages of As-Sawiya and al-Lubban al-Sharqiya, near Nablus, the official Palestinian news agency WAFA reported. According to local sources, Israeli forces threatened students and teachers […]
What is the Lions' Den and Why Now: 'The Palestinian View' – with Ramzy Baroud (VIDEO)
In his latest episode of 'The Palestinian View', Ramzy Baroud discusses the political context that led to the formation of the Lions' Den group in Nablus, leading to a growing armed rebellion throughout the West […]
Israel to Annex Palestinian Land for Settlement Expansion near Nablus
Israel will annex swathes of Palestinian land near Nablus in the occupied West Bank for the expansion of an illegal Israeli settlement, The New Arab reported. Ghassan Daghlas, who monitors settlement activities in the northern […]
Jewish Settlers Attack Palestinian Women Harvesting Olive Crops near Nablus
Jewish settlers on Tuesday attacked Palestinian women harvesting their olive crops in the village of Qaryout, near Nablus, the official Palestinian news agency WAFA reported. Ghassan Daghlas, a local official, said that settlers attempted to […]
Israeli Forces Kill Two Palestinian Members of Civil Defense near Nablus
Israeli occupation forces shot and killed two Palestinians on Thursday night near the Huwwara military checkpoint to the south of Nablus, the official Palestinian news agency WAFA reported. The Palestinian Ministry of Health said in […]
The Lions' Den Offers Details of What Has Actually Transpired in Nablus
On Monday, Israeli forces killed five Palestinians and wounded at least 21 in a major assault on the Old City of Nablus, in the occupied West Bank. Two completely different narratives emerged regarding the events […]


DuPage County Judge Kenneth Popejoy has been cited in a complaint filed Friday by the Illinois Judicial Inquiry Board that says he brought his "judicial office into disrepute" during a reckless-driving incident earlier this year.
The board is seeking a disciplinary hearing with the Illinois Courts Commission on the matter and is asking that action "deemed appropriate" be taken.
Popejoy was fined $500 and sentenced to six months of court supervision, a lesser form of probation, after pleading guilty to reckless driving on Aug. 3.
He originally was charged with failure to reduce speed to avoid an accident and failure to give information after striking an unattended motor vehicle on June 29.
According to court documents, Popejoy struck a parked vehicle in Glen Ellyn, then drove from the scene while the passenger-side front tire of his 2003 Jeep was nearly off the rim. En route to his home in Wheaton, he disobeyed two stop signs and caused a 13-year-old to jump out of the way to avoid being struck.
Popejoy, 59, who serves in the chancery division, has been a circuit judge since 2004 after serving as an associate judge since 1997.
Popejoy has declined to speak publicly on the case, citing a possible state hearing. After he was sentenced, Popejoy issued a statement saying he acted "in a moment of panic" and his decision to leave the scene and drive home "was an error in judgment that was inconsistent with my values and how I have lived my life."


This heavy metal band orientated in the sleepy Berkshire town of Sandhurst. So not the sort of place you'd expect to give birth to a heavy and thrash metal band! Primitai have combined influences from the 80s Glam Metal revival and classic hard rock and metal bands, such as Iron Maiden, AC/DC and Motley Crue. Creating melodic, rhythmic and powerful heavy metal. Whilst giving it a unique, modern twist to give them a unique edge in the resurgence of 80s Glam and Heavy Metal.
I was lucky to get two for one tonight, with two amazing heavy metal bands playing at the Black Heart, which is one of my favourite venues in Camden. Both fall into the broad heavy, power and thrash metal genres, but each take their own influences, styles and give it a unique twist.
Primitai are a treat to watch, their energy is undeniable and instantly infectious, spreading through that dark crowded room. With their electric fast riffs pumping out into the night, creating to true thrash metal experience. I simply adored their whole style.
They're high energy, hard-core, and die hard metal-heads. With a heavy yet fast beat, getting the crowd moving. Then the heavy crash of the drums layered with the psychedelic riffs created a mind-blowing experience.


MATH GETS REAL Chaos theory explains why it's so difficult to predict weather — a small change in conditions at any point in time can have a large effect on future conditions in, for example, a hurricane's trajectory (Tropical Storm Harvey shown).
Weird Math sets out to "reveal the strange connections between math and everyday life." The book fulfills that laudable goal, in part. At times, teenage math prodigy Agnijo Banerjee and his tutor, science writer David Darling, find ways to make complex math relatable, like linking chaos theory to weather forecasting and virtual reality. But there's a tension between precision and accessibility, and the authors don't always find the sweet spot.


THREE TIFFANY CUT DOUBLE OLD FASHIONED GLASSES. THIS IS A GORGEOUS SET OF THREE DOUBLE OLD FASHIONED GLASSES MADE BY TIFFANY AND COMPANY!!
THESE WONDERFUL GLASSES ARE VERY CLASSY AND HAVE HORIZONTAL AND VERTICAL CUT DESIGNS IN A TARTAN PLAID PATTERN. THEY MEASURE 3 5/8TH " HIGH AND 3 1/4" IN DIAMETER. THEY ARE IN EXCELLENT CONDITION WITH NO CHIPS OR DAMAGE OF ANY KIND.
ENTERTAIN YOUR GUEST WITH STYLE THIS HOLIDAY SEASON!! The item "TIFFANY,' TARTAN PLAID' CUT CRYSTAL DOUBLE OLD FASHIONED GLASSES (3) GORGEOUS" is in sale since Friday, November 2, 2018. This item is in the category "Pottery & Glass\Glass\Glassware\Cut Glass\Contemporary".
The seller is "maandpatina" and is located in Thomasville, Georgia. This item can be shipped worldwide.


The above words summarize in a line what traveling fundamentally is. But to live those words to the fullest means taking a leap of faith into the unknown. That unknown, which unfolds so many experiences that transform us and eventually become LIFE itself. On a lighter note, cutting off yourself from the daily chores, slowing down a bit and surrounding yourself with new people and being at new places can be an engaging and fulfilling experience. There are many reasons one can find to travel, and below are my 9 reasons to travel.
Have u seen the curiosity in the eyes of a toddler? The enthusiasm to learn something new and the happiness derived from it, those "first-time" experiences opening up to life's enrichment and the sense of contentment. Well, that's just the beginning. When you travel aimlessly without any preconceived notions about the places or the people, you are fueling that child within you and that is rewarding by itself.
We humans tend to live life in more or less a scheduled timetable to perform duties which makes our growth a limited function of time. Nothing can be as good as Traveling to know yourself better as you get lost into the unknown from the boundaries of your daily life. When you break the shackles of that confinement you live in, you get to see yourself evolved. Traveling gives that opportunity, time and situation which opens you up to new challenges and experiences which can help you identify your hidden strengths and helps you overcome some innate fears. Traveling also provides the platform to contemplate which all of us need at some point in life.
Yes that's right. Meeting new people, experiencing new culture, new surroundings and circumstances, enjoying new cuisines and visiting new places. All these add up to your richness. The friends you would make around the world, the experienced gourmet that you would turn into, the intellectually affluent person you would be and the fearless explorer you would become, all this would add up strongly in your contentment kitty. And the story telling that follows these experiences would even put to shame a man with really deep pockets who can't travel.
Most of the times we are confined within the boundaries drawn by our mind's reach which is in turn moderated by our perceptions and limited experiences as compared to the diversity globally. When you travel, you come out of those perceptions and get exposed to new people and their mind set which can add new dimensions to your thinking. When you meet new people and share your thought process on varied topics, you may get to know so much which your mind had never perceived to be true. This can be rewarding as it enriches your thoughts and makes you broad minded. And the view of life from a totally new pair of eyes, where you are an outsider, can be highly rewarding.
Traveling surely takes you out of your comfort zone. The impromptu changes in itinerary, new place, new people, sometimes new language to deal with, all this and much more gives you the ability to slow you down to a pace where you start identifying your traits in those situations. Every new experience makes you more stronger and ready to take on any challenges. Be it navigating for a good meal in France, or a walk in the most notorious street in Katmandu, to even getting invited for a conversation over Chai in Grand Bazaar in Istanbul, it all grooms you to be a one man army with unshakeable confidence.
Sometimes life becomes such a rut that you need a break to come out of it with zest. Traveling to places which are rich in its history, heritage, culture, food can be awe inspiring. A visit to a European museum, to see the hieroglyphs depicting the life of a Pharaoh, to see a war memorial or may be even a walk in the rain in a romantic city, can be sometimes be a medium to inspire you to do things that you always dreamt off.
Can be philosophical, but it's a fact that many people even in their mid age are not able to identify what makes them happy and what they really want to do in life. Also, most of the people can't cross over the social stigma to live life as they would have dreamt off, but traveling can sometimes act as the catalyst in the process to introspect on the strengths and weaknesses of a human and to eventually identify something that contents you.
With the fast paced life we live in, everyone burns out at a point in time and we all need a relaxing and pampering break. Travel can surely give you that much needed time to reorganize. Be it a trek in the Himalayan ranges, or a country side stay in some of the most beautiful villages of France, travel can surely help you unwind.
Have you been an audience to someone's storytelling of an adventure trip filled with bungee jumping or sky diving or trek in the wilderness, or some explorative travel like a journey on the Silk Road or that of a visit to Morocco, how about actually living one of those tales? Usually it's the human nature to seek the unknown and being experimentative. Travel can surely fuel your adventurous and explorative soul.
So what are you waiting for, book a ticket, go somewhere you haven't gone, meet new people, explore new places, live life and tell a tale worth an adventure.
How very well said! I came back to read your blog after a while. Agree every bit of what you have written about travelling. Sometimes, you are so addicted to travelling that you want to travel with any little break you get! New perspective, new experiences, the way you see other people and cultures changes! I would love to read about your experiences of the Everest Base Camp! Is it in the pipeline?
I think that travelers definitely become more independent as they take long term trips. You have to because if not then you'll have no one to take care of yourself.